Coverage for src/fast_minimum_variance/minvar_problem.py: 100%
228 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-04 15:49 +0000
1"""Minimum-variance solver: primal asset elimination with dual-feasibility check."""
3from collections.abc import Callable
4from dataclasses import dataclass
5from typing import Any
7import numpy as np
8from cvx.linalg import DenseOperator, FactorOperator, GramOperator, SumOperator
9from scipy.linalg import solve as spd_solve
10from scipy.sparse.linalg import LinearOperator, cg
12from ._base import _BaseProblem
15@dataclass(frozen=True)
16class _MinVarProblem(_BaseProblem):
17 """Minimum-variance portfolio solver via primal-dual active-set iteration.
19 Solves::
21 min (1-alpha)||X w||^2 + alpha*(||X||_F^2/N)*||w||^2 - rho*mu^T w
22 s.t. B w = c, w >= 0
24 where ``B`` is a ``(p, N)`` balance system with full row rank on every
25 active set the loop visits. The default (``B = None``) is the budget
26 constraint ``1^T w = 1`` (``p = 1``).
28 Each inner step solves the equality-constrained subproblem over the current
29 active asset set. Stationarity gives ``2*Sigma_a*w_a = B_a^T lambda + rho*mu_a``
30 where ``Sigma_a = (1-alpha)*X_a^T X_a + ridge*I``. Solving the ``n_a x n_a``
31 SPD system ``Sigma_a V = B_a^T`` (``p`` right-hand sides, plus
32 ``Sigma_a v_mu = mu_a`` when ``rho != 0``) and recovering ``lambda`` from the
33 ``p x p`` Schur system ``(B_a V) lambda = c`` avoids the indefinite
34 ``(n_a+p) x (n_a+p)`` saddle-point system entirely. The outer primal-dual
35 loop enforces ``w >= 0`` and terminates when both primal and dual feasibility
36 hold simultaneously.
38 Use ``alpha = N/(N+T)`` for Ledoit-Wolf shrinkage intensity::
40 T, N = X.shape
41 w, iters = Problem(X, alpha=N/(N+T)).solve_kkt()
43 Examples:
44 >>> import numpy as np
45 >>> from fast_minimum_variance import Problem
46 >>> X = np.random.default_rng(0).standard_normal((100, 5))
47 >>> w, iters = Problem(X).solve_kkt()
48 >>> float(round(w.sum(), 6))
49 1.0
50 >>> bool((w >= 0).all())
51 True
52 """
54 B: np.ndarray | None = None # (p, N) balance system: B w = c; None = budget 1^T w = 1
55 c: np.ndarray | None = None # (p,) balance targets
57 def __post_init__(self) -> None:
58 """Validate balance-system shapes on top of the base checks."""
59 super().__post_init__()
60 if (self.B is None) != (self.c is None):
61 raise ValueError("B and c must be supplied together") # noqa: TRY003
62 if self.B is not None:
63 c = self.c
64 assert c is not None # noqa: S101
65 if self.B.ndim != 2 or self.B.shape[1] != self.n:
66 raise ValueError(f"B must have shape (p, {self.n}), got {self.B.shape}") # noqa: TRY003
67 if c.shape != (self.B.shape[0],):
68 raise ValueError(f"c must have shape ({self.B.shape[0]},), got {c.shape}") # noqa: TRY003
70 # ------------------------------------------------------------------
71 # Balance-system helpers (budget is the p = 1 special case)
72 # ------------------------------------------------------------------
74 @property
75 def _p(self) -> int:
76 """Number of balance constraints (1 for the default budget)."""
77 return 1 if self.B is None else int(self.B.shape[0])
79 def _c_vec(self) -> np.ndarray:
80 """Balance right-hand side ``c`` (the budget gives ``[1.0]``)."""
81 return np.ones(1) if self.c is None else self.c
83 def _balance_rows(self, active: np.ndarray) -> np.ndarray:
84 """Return ``B_a``, the balance system restricted to active assets, shape ``(p, n_a)``."""
85 if self.B is None:
86 return np.ones((1, int(active.sum())))
87 return self.B[:, active]
89 def _recover_balance(self, v_eq: np.ndarray, v_mu: np.ndarray | None, b_a: np.ndarray) -> np.ndarray:
90 """Recover ``w_a`` from the Schur reduction of the balance system.
92 Given ``V = Sigma_a^{-1} B_a^T`` (columns of ``v_eq``) and optionally
93 ``v_mu = Sigma_a^{-1} mu_a``, stationarity ``2*Sigma_a*w = B_a^T lambda
94 + rho*mu_a`` and feasibility ``B_a w = c`` pin the multiplier through
95 the ``p x p`` SPD Schur system ``(B_a V) eta = c - rho/2 * B_a v_mu``
96 with ``eta = lambda/2``, so ``w = V eta + rho/2 * v_mu``.
97 """
98 schur = b_a @ v_eq # (p, p)
99 rhs = self._c_vec().astype(np.float64)
100 if v_mu is not None:
101 rhs = rhs - 0.5 * self.rho * (b_a @ v_mu)
102 eta = np.linalg.solve(schur, rhs) if self._p > 1 else rhs / schur[0, 0]
103 w = v_eq @ eta
104 if v_mu is not None:
105 w = w + 0.5 * self.rho * v_mu
106 return w
108 # ------------------------------------------------------------------
109 # Shared helpers used by both active-set loop variants
110 # ------------------------------------------------------------------
112 def _compute_gradient(self, w: np.ndarray) -> np.ndarray:
113 """Return the full objective gradient at w, including rho*mu adjustment."""
114 data_grad = (self.X.T @ (self.X @ w)) / self.t
115 if self.target_lr is not None:
116 bar_lam, U_k, delta_k = self.target_lr # noqa: N806
117 tgt_w = bar_lam * w + U_k @ (delta_k * (U_k.T @ w))
118 grad = 2.0 * ((1 - self.alpha) * data_grad + self.alpha * tgt_w)
119 elif self.target is not None:
120 grad = 2.0 * ((1 - self.alpha) * data_grad + self.alpha * self.target @ w)
121 else:
122 grad = 2.0 * data_grad
123 if self.rho != 0.0 and self.mu is not None:
124 grad = grad - self.rho * self.mu
125 result: np.ndarray = grad
126 return result
128 @staticmethod
129 def _primal_drop(w_a: np.ndarray, asset_active: np.ndarray, tol: float) -> bool:
130 """Drop negative-weight assets from active set in-place; return True if any dropped."""
131 if not np.any(w_a < -tol):
132 return False
133 idx = np.where(asset_active)[0]
134 strong = w_a < -10 * tol
135 if np.any(strong):
136 asset_active[idx[strong]] = False
137 else:
138 asset_active[idx[np.argmin(w_a)]] = False
139 return True
141 def _dual_add(self, grad: np.ndarray, asset_active: np.ndarray, tol: float) -> int:
142 """Return index of excluded asset that violates KKT dual condition, or -1 if none.
144 The multiplier is estimated from the active gradient: for the budget the
145 stationary ``lambda`` is a location estimate of ``g_a`` (median for
146 robustness on larger sets); for a general balance system it is the
147 least-squares solution of ``B_a^T lambda = g_a``. The bound multiplier
148 estimate is then ``nu = grad - B^T lambda``, which must be non-negative
149 on excluded assets at the optimum.
150 """
151 excluded = ~asset_active
152 if not excluded.any():
153 return -1
154 g_a = grad[asset_active]
155 if self.B is None:
156 lambda_ = np.median(g_a) if g_a.size > 5 else g_a.mean()
157 nu = grad - lambda_
158 else:
159 b_a = self.B[:, asset_active]
160 lam, *_ = np.linalg.lstsq(b_a.T, g_a, rcond=None)
161 nu = grad - self.B.T @ lam
162 idx_ex = np.where(excluded)[0]
163 j = idx_ex[np.argmin(nu[excluded])]
164 return int(j) if nu[j] < -tol else -1
166 # ------------------------------------------------------------------
167 # Outer loop: primal elimination + dual feasibility check
168 # ------------------------------------------------------------------
169 def _constraint_active_set(
170 self,
171 solve_fn: Callable[[np.ndarray], tuple[np.ndarray, int]],
172 tol: float = 1e-6,
173 max_iter: int = 10_000,
174 ) -> tuple[np.ndarray, int, int]:
175 """Run the primal-dual active-set loop enforcing ``w >= 0``.
177 Calls ``solve_fn(active_mask)`` repeatedly. The *primal step* drops assets
178 with negative weights; the *dual step* re-adds any excluded asset whose KKT
179 gradient condition is violated. Terminates when both conditions hold
180 simultaneously, which together with stationarity is sufficient for global
181 optimality.
182 """
183 n = self.n
184 asset_active = np.ones(n, dtype=bool)
185 total_inner_iters = 0
186 outer_steps = 0
187 prev_active = None
188 w = np.zeros(n)
190 for _ in range(max_iter):
191 if prev_active is not None and np.array_equal(prev_active, asset_active):
192 break # pragma: no cover - structurally unreachable safety guard
193 prev_active = asset_active.copy()
195 w_a, step_iters = solve_fn(asset_active)
196 outer_steps += 1
197 total_inner_iters += step_iters
199 if self._primal_drop(w_a, asset_active, tol):
200 continue
202 w = np.zeros(n)
203 w[asset_active] = w_a
205 j = self._dual_add(self._compute_gradient(w), asset_active, tol)
206 if j < 0:
207 break
208 asset_active[j] = True
210 return w, outer_steps, total_inner_iters
212 # ------------------------------------------------------------------
213 # Inner steps
214 # ------------------------------------------------------------------
216 def _kkt_step(self, active: np.ndarray, x0: np.ndarray | None = None) -> tuple[np.ndarray, int]: # noqa: ARG002
217 """Solve the reduced SPD system directly; return ``(w_a, 1)``.
219 Stationarity gives ``2*Sigma_a*w_a = B_a^T lambda + rho*mu_a``. A single
220 solve with ``p`` RHS columns yields ``V = Sigma_a^{-1} B_a^T`` (plus
221 ``v_mu = Sigma_a^{-1} mu_a`` when ``rho != 0``); the balance system then
222 pins ``lambda`` through the ``p x p`` Schur solve of
223 :meth:`_recover_balance` (for the budget this reduces to
224 ``lambda = 2*(1 - rho/2 * sum(v_mu)) / sum(v1)``).
226 When ``alpha=1`` and ``target_lr`` is set the system is purely the RMT
227 target ``T0 = bar_lam*I + U_k diag(delta_k) U_k^T``. The Woodbury
228 identity gives the exact inverse in O(n_a*k + k^3) without CG iterations:
229 ``T0^{-1} b = b/bar_lam - U_k_a W^{-1}(U_k_a^T b)/bar_lam^2``
230 where ``W = diag(1/delta_k) + U_k_a^T U_k_a / bar_lam``.
231 """
232 b_a = self._balance_rows(active)
233 tilt = self.rho != 0.0 and self.mu is not None
234 rhs = np.column_stack([b_a.T, self.mu[active]]) if tilt and self.mu is not None else b_a.T
236 # Woodbury direct solve: O(n_a*k + k^3) for alpha=1, RMT target. The target
237 # T0 = bar_lam*I + U_k diag(delta_k) U_k^T is a diagonal-plus-low-rank operator,
238 # so its active-block inverse is exactly cvx-linalg's FactorOperator.solve_free.
239 if self.alpha == 1.0 and self.target_lr is not None:
240 bar_lam, U_k, delta_k = self.target_lr # noqa: N806
241 t0 = FactorOperator(np.full(U_k.shape[0], bar_lam), U_k, np.diag(delta_k))
242 sols = np.asarray(t0.solve_free(np.flatnonzero(active), rhs))
243 else:
244 x_a = self.X[:, active]
245 if self.target is None:
246 sigma = (x_a.T @ x_a) / self.t
247 else:
248 sigma = (1.0 - self.alpha) * (x_a.T @ x_a) / self.t + self.alpha * self.target[np.ix_(active, active)]
249 sols = spd_solve(sigma, rhs, assume_a="pos")
251 v_mu = sols[:, -1] if tilt else None
252 return self._recover_balance(sols[:, : self._p], v_mu, b_a), 1
254 def _cvxpy_constraints(self, w: Any, cp: Any) -> list[Any]:
255 """Return balance-equality and long-only inequality constraints for CVXPY."""
256 if self.B is not None:
257 return [self.B @ w == self.c, w >= 0]
258 return [cp.sum(w) == 1, w >= 0]
260 def _system_operator(self) -> SumOperator:
261 """Build ``Sigma = (1-alpha)/T * X^T X + alpha * T0`` as a cvx-linalg operator.
263 A :class:`~cvx.linalg.SumOperator` of the data Gram term and, when present,
264 the target term (a :class:`~cvx.linalg.FactorOperator` for a low-rank RMT
265 target, else a :class:`~cvx.linalg.DenseOperator`). The full-universe
266 operators are sliced to the active set via ``apply_free``; nothing is
267 formed at ``n x n``. Without a target the data term carries the full weight.
268 """
269 has_target = self.target_lr is not None or self.target is not None
270 c_data = (1.0 - self.alpha) if has_target else 1.0
271 terms: list[tuple[float, Any]] = [(c_data / self.t, GramOperator(self.X))]
272 if self.target_lr is not None:
273 bar_lam, u_k, delta_k = self.target_lr
274 terms.append((self.alpha, FactorOperator(np.full(u_k.shape[0], bar_lam), u_k, np.diag(delta_k))))
275 elif self.target is not None:
276 terms.append((self.alpha, DenseOperator(self.target)))
277 return SumOperator(terms)
279 @staticmethod
280 def _free_matvec(sigma: SumOperator, active_idx: np.ndarray) -> Callable[[np.ndarray], np.ndarray]:
281 """Return the free-block action ``v -> Sigma[A, A] v`` with the slice hoisted out.
283 When cvx-linalg exposes ``restricted`` (>= 0.9.6), the pre-sliced free-block
284 operator is built once here and its plain ``matvec`` is returned. Calling
285 ``apply_free(idx, v)`` per CG iteration instead re-gathers the operator's
286 storage (the Gram factor columns) on every call, which costs an order of
287 magnitude more wall clock at identical iteration counts. The fallback keeps
288 older cvx-linalg releases working.
289 """
290 restricted = getattr(sigma, "restricted", None)
291 if restricted is not None:
292 try:
293 restricted_op = restricted(active_idx)
294 except NotImplementedError:
295 restricted_op = None
296 if restricted_op is not None:
297 matvec: Callable[[np.ndarray], np.ndarray] = restricted_op.matvec
298 return matvec
299 return lambda v: sigma.apply_free(active_idx, v)
301 def _cg_step(self, active: np.ndarray, x0: np.ndarray | None = None) -> tuple[np.ndarray, int]:
302 """Solve the reduced SPD system via matrix-free CG; return ``(w_a, iters)``.
304 Runs conjugate gradients over the active-set system operator
305 (:meth:`_system_operator`), restricted to the active set once per step so
306 the reduced matvec is ``O(n_a T)`` rather than ``O(n T)``, without ever
307 forming ``Sigma_a`` explicitly. Low-rank and dense targets share one path.
309 Args:
310 active: Boolean mask selecting the active asset subset.
311 x0: Optional initial guess for the first CG solve (warm start).
312 When provided it must have shape ``(active.sum(),)``.
313 """
314 n_a = int(active.sum())
315 sigma = self._system_operator()
316 active_idx = np.flatnonzero(active)
317 free_matvec = self._free_matvec(sigma, active_idx)
318 count = [0]
320 def matvec(v: np.ndarray) -> np.ndarray:
321 """Apply Sigma_a to v via the pre-sliced free-block operator."""
322 count[0] += 1
323 return free_matvec(v)
325 op = LinearOperator((n_a, n_a), matvec=matvec, dtype=np.float64) # ty:ignore[missing-argument, parameter-already-assigned, unknown-argument]
327 b_a = self._balance_rows(active)
328 # x0 approximates the final w, which is proportional to the single
329 # solve column only in the budget case; skip the guess for p > 1.
330 guess = x0 if self._p == 1 else None
331 v_eq = np.column_stack([cg(op, b_a[j], x0=guess)[0] for j in range(self._p)])
332 v_mu = cg(op, self.mu[active], x0=guess)[0] if self.rho != 0.0 and self.mu is not None else None
333 return self._recover_balance(v_eq, v_mu, b_a), count[0]
335 def _pcg_step(self, active: np.ndarray, x0: np.ndarray | None = None) -> tuple[np.ndarray, int]:
336 """Solve the reduced SPD system via PCG with RMT preconditioner; return (w_a, iters).
338 The system matrix is the oracle-LW covariance (using self.alpha and self.target).
339 The preconditioner P = T0^RMT is applied via the Woodbury identity:
340 P^{-1} v = (1/bar_lam) v + U_k diag(1/lambda_k - 1/bar_lam) U_k^T v
341 costing O(n_a * k) per application. Requires self.pcg_lr to be set.
342 """
343 n_a = int(active.sum())
345 # System matvec — the same active-set operator as _cg_step, sliced once.
346 sigma = self._system_operator()
347 active_idx = np.flatnonzero(active)
348 free_matvec = self._free_matvec(sigma, active_idx)
349 count = [0]
351 def matvec(v: np.ndarray) -> np.ndarray:
352 """Apply the active-set system matrix Sigma_a to v."""
353 count[0] += 1
354 return free_matvec(v)
356 op = LinearOperator((n_a, n_a), matvec=matvec, dtype=np.float64) # ty:ignore[missing-argument, parameter-already-assigned, unknown-argument]
358 # Preconditioner P^{-1}: Woodbury inverse of T0^RMT restricted to active set
359 pcg_lr = self.pcg_lr
360 if pcg_lr is None: # pragma: no cover - defensive; solve_pcg validates pcg_lr upfront (see _base.py)
361 raise RuntimeError("_pcg_step called without pcg_lr") # noqa: TRY003
362 bar_lam_p, U_k_p, delta_k_p = pcg_lr # noqa: N806
363 U_k_a_p = U_k_p[active, :] # noqa: N806 # (n_a, k)
364 inv_coeff = 1.0 / (bar_lam_p + delta_k_p) - 1.0 / bar_lam_p # (k,) negative
366 def precond(v: np.ndarray) -> np.ndarray:
367 """Apply P^{-1} to v via the Woodbury identity."""
368 result: np.ndarray = (1.0 / bar_lam_p) * v + U_k_a_p @ (inv_coeff * (U_k_a_p.T @ v))
369 return result
371 M_op = LinearOperator((n_a, n_a), matvec=precond, dtype=np.float64) # ty:ignore[missing-argument, parameter-already-assigned, unknown-argument] # noqa: N806
373 b_a = self._balance_rows(active)
374 guess = x0 if self._p == 1 else None
375 v_eq = np.column_stack([cg(op, b_a[j], x0=guess, M=M_op)[0] for j in range(self._p)])
376 return self._recover_balance(v_eq, None, b_a), count[0]
378 def _constraint_active_set_warm(
379 self,
380 solve_fn: Callable[..., tuple[np.ndarray, int]] | None = None,
381 tol: float = 1e-6,
382 max_iter: int = 10_000,
383 warm_start: tuple[np.ndarray, np.ndarray] | None = None,
384 ) -> tuple[np.ndarray, int, int, np.ndarray, np.ndarray | None]:
385 """Active-set loop with warm-starting; returns ``(w, iters, active, w_full)``.
387 Generalises ``_constraint_active_set``: accepts an initial active set
388 and passes the previous iterate as a starting guess to the inner solver.
389 Solvers that support an initial guess (CG via ``_cg_step``) benefit from
390 both the warm active set and the x0; direct solvers (KKT via
391 ``_kkt_step``) accept and silently ignore the x0, profiting only from
392 the warm active set.
394 Args:
395 solve_fn: Inner solver callable ``(active, x0=None) -> (w_a, iters)``.
396 Defaults to ``self._cg_step``.
397 tol: Primal feasibility tolerance; assets with weight below ``-tol``
398 are dropped from the active set.
399 max_iter: Maximum number of outer active-set iterations.
400 warm_start: Optional ``(active_mask, w_full)`` from a previous call.
401 ``active_mask`` is a boolean array of length ``n``;
402 ``w_full`` is the full ``n``-vector of weights.
404 Returns:
405 ``(w, total_iters, final_active, w_full)`` — solution, cumulative
406 iteration count, final active-set mask, and full weight vector
407 suitable as the ``warm_start`` argument for the next solve.
408 """
409 if solve_fn is None:
410 solve_fn = self._cg_step
411 n = self.n
413 if warm_start is not None:
414 asset_active, last_w_full = warm_start
415 asset_active = asset_active.copy()
416 else:
417 asset_active = np.ones(n, dtype=bool)
418 last_w_full = None
420 total_inner_iters = 0
421 outer_steps = 0
422 prev_active = None
423 w = np.zeros(n)
425 for _ in range(max_iter):
426 if prev_active is not None and np.array_equal(prev_active, asset_active):
427 break # pragma: no cover - structurally unreachable safety guard
428 prev_active = asset_active.copy()
430 x0 = None
431 if last_w_full is not None:
432 sub = last_w_full[asset_active]
433 s = sub.sum()
434 if s > 1e-12:
435 x0 = sub / s
437 w_a, step_iters = solve_fn(asset_active, x0=x0)
438 outer_steps += 1
439 total_inner_iters += step_iters
441 if self._primal_drop(w_a, asset_active, tol):
442 continue
444 w = np.zeros(n)
445 w[asset_active] = w_a
446 last_w_full = w.copy()
448 j = self._dual_add(self._compute_gradient(w), asset_active, tol)
449 if j < 0:
450 break
451 asset_active[j] = True
453 return w, outer_steps, total_inner_iters, asset_active.copy(), last_w_full
455 def solve_cg_warm(
456 self,
457 *,
458 project: bool = True,
459 warm_start: tuple[np.ndarray, np.ndarray] | None = None,
460 ) -> tuple[np.ndarray, int, int, tuple[np.ndarray, np.ndarray | None]]:
461 """Solve via matrix-free CG with warm-starting.
463 Like ``solve_cg`` but accepts and returns warm-start state so that a
464 sequence of related problems (e.g. an efficient-frontier sweep over
465 many ``rho`` values) can chain solves together. Adjacent problems share
466 a similar active set and similar solution, so subsequent solves need
467 far fewer outer iterations and CG steps.
469 Args:
470 project: Clip weights to ``[0, ∞)`` and renormalize to sum to 1.
471 warm_start: ``(active_mask, w_full)`` returned by a previous call,
472 or ``None`` for a cold start.
474 Returns:
475 ``(w, outer_steps, inner_iters, warm_state)`` — weight vector,
476 number of outer active-set steps, cumulative CG iteration count,
477 and warm state for the next call in the sequence.
479 Examples:
480 >>> import numpy as np
481 >>> from fast_minimum_variance import Problem
482 >>> X = np.random.default_rng(0).standard_normal((100, 5))
483 >>> mu = np.ones(5) * 0.01
484 >>> warm = None
485 >>> for rho in [0.0, 0.5, 1.0]:
486 ... p = Problem(X, rho=rho, mu=mu)
487 ... w, outer, inner, warm = p.solve_cg_warm(warm_start=warm)
488 >>> float(round(w.sum(), 10))
489 1.0
490 >>> bool((w >= 0).all())
491 True
492 """
493 w, outer, inner, final_active, final_w = self._constraint_active_set_warm(
494 solve_fn=self._cg_step, warm_start=warm_start
495 )
496 if project:
497 w = self._clip_and_renormalize(w)
498 return w, outer, inner, (final_active, final_w)
500 def solve_kkt_warm(
501 self,
502 *,
503 project: bool = True,
504 warm_start: tuple[np.ndarray, np.ndarray] | None = None,
505 ) -> tuple[np.ndarray, int, tuple[np.ndarray, np.ndarray | None]]:
506 """Solve via direct KKT factorisation with active-set warm-starting.
508 Like ``solve_kkt`` but accepts and returns warm-start state for chaining
509 a sequence of related problems. KKT is a direct solver, so only the
510 active-set mask is warm-started (not an iterative initial guess); the
511 benefit is fewer outer primal-dual iterations when consecutive problems
512 share a similar active set.
514 Args:
515 project: Clip weights to ``[0, ∞)`` and renormalize to sum to 1.
516 warm_start: ``(active_mask, w_full)`` returned by a previous call,
517 or ``None`` for a cold start.
519 Returns:
520 ``(w, outer_steps, warm_state)`` — weight vector, number of outer
521 active-set steps, and warm state for the next call in the sequence.
523 Examples:
524 >>> import numpy as np
525 >>> from fast_minimum_variance import Problem
526 >>> X = np.random.default_rng(0).standard_normal((100, 5))
527 >>> mu = np.ones(5) * 0.01
528 >>> warm = None
529 >>> for rho in [0.0, 0.5, 1.0]:
530 ... p = Problem(X, rho=rho, mu=mu)
531 ... w, outer, warm = p.solve_kkt_warm(warm_start=warm)
532 >>> float(round(w.sum(), 10))
533 1.0
534 >>> bool((w >= 0).all())
535 True
536 """
537 w, outer, _inner, final_active, final_w = self._constraint_active_set_warm(
538 solve_fn=self._kkt_step, warm_start=warm_start
539 )
540 if project:
541 w = self._clip_and_renormalize(w)
542 return w, outer, (final_active, final_w)
544 # ------------------------------------------------------------------
545 # Budget-specific overrides
546 # ------------------------------------------------------------------
548 def _clip_and_renormalize(self, w: np.ndarray) -> np.ndarray: # type: ignore[override]
549 """Project onto the budget simplex; identity when a balance system is set.
551 Renormalising by the weight sum would break a general ``B w = c``, and
552 the active-set loop already exits primal-feasible, so balance-system
553 solves return the iterate unchanged.
554 """
555 if self.B is not None:
556 return w
557 return _BaseProblem._clip_and_renormalize(w)