Coverage for src/cvxcla/cla.py: 100%
176 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-11 09:58 +0000
1"""Markowitz implementation of the Critical Line Algorithm.
3This module provides the CLA class, which implements the Critical Line Algorithm
4as described by Harry Markowitz and colleagues. The algorithm computes the entire
5efficient frontier by finding all turning points, which are the points where the
6set of assets at their bounds changes.
7"""
9import logging
10from dataclasses import dataclass, field
11from functools import cached_property
12from typing import NamedTuple
14import numpy as np
15from numpy.typing import NDArray
17from ._events import event_ratios, ineq_event_ratios
18from ._kkt import active_set, solve_kkt
19from ._projection import project_feasible
20from .first import first_vertex_lp, init_algo
21from .operators import DenseCovariance, QuadraticForm
22from .pathtracer import InequalityConstrained, trace
23from .types import Frontier, FrontierPoint, TurningPoint
25# A genuinely rank-deficient free block has a reciprocal condition number at
26# round-off level (~1e-16); a well-posed or merely near-degenerate block sits
27# many orders above it (>= ~1e-4 across the degeneracy sweep in
28# experiments/degeneracy_boundary.py). The 1e-12 cut sits in the wide gap between
29# the two and is the conventional numerical-singularity scale.
30_RCOND_FLOOR = 1e-12 # pragma: no mutate
33class _Segment(NamedTuple):
34 """The affine critical-line segment valid at one turning point.
36 Bundles the affine path ``w(lam) = r_alpha + lam * r_beta``, the multiplier
37 gradients ``gamma``/``delta`` that drive the leave-a-bound events, and the
38 active-set masks the event scan needs. This is what ``CLA.segment`` returns
39 to the generic path tracer.
41 For general inequality constraints ``G w <= h`` the segment also carries the
42 affine inequality multipliers ``eta_alpha + lam * eta_beta`` (one entry per
43 inequality row; meaningful for *active* rows, which release when the
44 multiplier crosses zero) and the active-row mask ``active_ineq``. The slacks
45 that drive an *inactive* row becoming active are recomputed from
46 ``r_alpha``/``r_beta`` directly in :func:`cvxcla._events.ineq_event_ratios`.
47 """
49 r_alpha: NDArray[np.float64]
50 r_beta: NDArray[np.float64]
51 gamma: NDArray[np.float64]
52 delta: NDArray[np.float64]
53 at_upper: NDArray[np.bool_]
54 at_lower: NDArray[np.bool_]
55 free_in: NDArray[np.bool_]
56 active_ineq: NDArray[np.bool_]
57 eta_alpha: NDArray[np.float64]
58 eta_beta: NDArray[np.float64]
61@dataclass(frozen=True)
62class CLA(InequalityConstrained):
63 """Critical Line Algorithm implementation based on Markowitz's approach.
65 This class implements the Critical Line Algorithm as described by Harry Markowitz
66 and colleagues. It computes the entire efficient frontier by finding all turning
67 points, which are the points where the set of assets at their bounds changes.
69 The algorithm starts with the first turning point (the portfolio with the highest
70 expected return) and then iteratively computes the next turning point with a lower
71 expected return until it reaches the minimum variance portfolio.
73 Attributes:
74 mean: Vector of expected returns for each asset.
75 covariance: Covariance matrix of asset returns, either as a plain
76 ``numpy`` array or as a ``CovarianceOperator`` backend
77 (see ``cvxcla.operators``).
78 lower_bounds: Vector of lower bounds for asset weights.
79 upper_bounds: Vector of upper bounds for asset weights.
80 a: Equality-constraint matrix ``A`` of ``A w = b`` (``m x n``). The
81 canonical case is the single all-ones budget row (``sum(w) = b``),
82 but an arbitrary equality system is supported: weighted single rows
83 and ``m > 1`` rows (e.g. budget plus sector- or factor-neutrality).
84 The all-ones budget (any right-hand side, including ``0`` for
85 dollar-neutral) uses the greedy first vertex of
86 :func:`cvxcla.first.init_algo`; a general ``A`` uses the
87 linear-programming first vertex of
88 :func:`cvxcla.first.first_vertex_lp`.
89 b: Equality-constraint right-hand side ``b`` (length ``m``); ``[1]`` for
90 the fully-invested budget, ``[0]`` for dollar-neutral, and so on.
91 g: Optional inequality-constraint matrix ``G`` of ``G w <= h``
92 (``p x n``), e.g. a group- or sector-exposure cap. ``None`` (the
93 default) means no inequality rows, recovering the equality-only
94 problem exactly. A ``>=`` constraint is expressed by negating both
95 ``g`` and ``h``. Each *active* row (held at equality ``g_i w = h_i``)
96 enters the reduced KKT system as an extra equality row, so the
97 covariance is still touched only through the ``QuadraticForm``
98 interface; box bounds remain a separate per-variable active set.
99 h: Optional inequality-constraint right-hand side ``h`` (length ``p``).
100 turning_points: List of turning points on the efficient frontier.
101 tol: Tolerance for numerical calculations.
102 logger: Logger instance for logging information and errors.
104 """
106 mean: NDArray[np.float64]
107 covariance: NDArray[np.float64] | QuadraticForm
108 lower_bounds: NDArray[np.float64]
109 upper_bounds: NDArray[np.float64]
110 a: NDArray[np.float64]
111 b: NDArray[np.float64]
112 g: NDArray[np.float64] | None = None
113 h: NDArray[np.float64] | None = None
114 turning_points: list[TurningPoint] = field(default_factory=list)
115 tol: float = 1e-5 # pragma: no mutate
116 logger: logging.Logger = field(default_factory=lambda: logging.getLogger(__name__))
118 @classmethod
119 def problem(cls, mean: NDArray[np.float64], covariance: NDArray[np.float64] | QuadraticForm) -> "ProblemBuilder":
120 """Start a fluent :class:`cvxcla.builder.ProblemBuilder` for this problem.
122 A readability convenience over the explicit constructor: chain
123 ``.long_only()``/``.budget()``/``.equality()``/``.inequality()`` and finish
124 with ``.trace()``. The builder maps one-to-one onto the constructor
125 arguments and adds no modelling power.
127 Args:
128 mean: Vector of expected returns of length ``n``.
129 covariance: Covariance matrix or ``QuadraticForm`` backend.
131 Returns:
132 A ``ProblemBuilder`` ready to accept constraints.
133 """
134 return ProblemBuilder(mean, covariance)
136 @cached_property
137 def covariance_operator(self) -> QuadraticForm:
138 """Return the covariance as a ``QuadraticForm`` backend.
140 A plain ``numpy`` covariance matrix is wrapped in ``DenseCovariance``;
141 an object already implementing the protocol is passed through. This is
142 the single point where the input form is normalised.
143 """
144 if isinstance(self.covariance, QuadraticForm):
145 return self.covariance
146 return DenseCovariance(self.covariance)
148 @property
149 def dimension(self) -> int:
150 """Number of assets ``n`` (the problem dimension for the path tracer)."""
151 return len(self.mean)
153 @cached_property
154 def _free_blocks_well_conditioned(self) -> bool:
155 """Whether every free-block solve along the trace is numerically safe.
157 Decided once, up front. By Cauchy's interlacing theorem every principal
158 submatrix of the symmetric PSD covariance is at least as well conditioned
159 as the whole matrix -- deleting rows/columns cannot decrease the smallest
160 eigenvalue nor increase the largest -- so the reciprocal condition number
161 of any free block is ``>=`` that of the full covariance. Hence if the full
162 covariance clears the singularity floor, no free block encountered along
163 the trace can be singular, and the per-turning-point conditioning guard in
164 :meth:`_emit` is provably never triggered. We then skip it, paying one
165 conditioning test here instead of one at every turning point (the latter
166 is a full eigendecomposition of the free block, as costly as the KKT solve,
167 so it otherwise dominates the trace). The up-front test computes the
168 reciprocal condition number of the full covariance once, via
169 :meth:`~cvx.linalg.SymmetricOperator.rcond_free`.
171 When the full covariance is itself near-singular (for example a sample
172 covariance from fewer observations than assets) this is ``False`` and the
173 per-step guard in :meth:`_emit` runs unchanged, preserving the degeneracy
174 diagnosis exactly.
175 """
176 full = np.arange(self.dimension)
177 return self.covariance_operator.rcond_free(full) >= _RCOND_FLOOR
179 def __post_init__(self) -> None:
180 """Initialize the CLA object and compute the efficient frontier.
182 This method is automatically called after initialization. It computes
183 the entire efficient frontier by finding all turning points, starting
184 from the first turning point (highest expected return) and iteratively
185 computing the next turning point with a lower expected return until
186 it reaches the minimum variance portfolio.
188 The actual walk is driven by the generic ``cvxcla.pathtracer.trace``
189 loop; this class supplies the portfolio-specific hooks (``begin``,
190 ``segment``, ``event_matrix``, ``step``, ``finish``) it calls.
192 The reduced KKT system at each turning point is solved by block
193 elimination: a single multi-RHS solve against the free covariance block
194 (via the covariance backend), covering the constraint columns and the
195 alpha and beta systems together so ``Sigma_FF`` is factorised once, then a
196 small Schur-complement solve ``A_F @ Sigma_FF^{-1} @ A_F.T`` over the
197 equality (and active inequality) rows. The covariance only enters through
198 the ``QuadraticForm`` interface, so structured backends (e.g.
199 ``FactorCovariance``) never materialise an n x n matrix.
201 Raises:
202 RuntimeError: If all variables are blocked, which would make the
203 system of equations singular.
204 ValueError: If the inequality matrix ``g`` and vector ``h`` have
205 mismatched or wrong shapes.
207 """
208 if self.g_matrix.shape[1] != self.dimension:
209 msg = f"g must have {self.dimension} columns, got shape {self.g_matrix.shape}"
210 raise ValueError(msg)
211 if self.h_vector.shape[0] != self.g_matrix.shape[0]:
212 msg = f"h must have {self.g_matrix.shape[0]} entries, got {self.h_vector.shape[0]}"
213 raise ValueError(msg)
214 trace(self)
216 def begin(self) -> tuple[float, TurningPoint]:
217 """Record the first turning point and start the trace at ``lambda = inf``.
219 Returns:
220 ``(inf, first_turning_point)``: the starting lambda bound and the
221 initial state for the path tracer.
222 """
223 first = self._first_turning_point()
224 self._append(first)
225 return np.inf, first
227 def segment(self, state: TurningPoint) -> _Segment:
228 """Solve the reduced KKT system for the critical-line segment at ``state``."""
229 at_upper, at_lower, free_in, fixed_weights = active_set(
230 state.free, state.weights, self.lower_bounds, self.upper_bounds, self.tol
231 )
232 r_alpha, r_beta, gamma, delta, eta_alpha, eta_beta = solve_kkt(
233 self.covariance_operator,
234 self.mean,
235 self.a,
236 self.b,
237 self.g_matrix,
238 self.h_vector,
239 free_in,
240 fixed_weights,
241 state.active_ineq,
242 )
243 return _Segment(
244 r_alpha, r_beta, gamma, delta, at_upper, at_lower, free_in, state.active_ineq, eta_alpha, eta_beta
245 )
247 def event_matrix(self, state: TurningPoint, segment: _Segment) -> NDArray[np.float64]: # noqa: ARG002
248 """Return the ``(n + p, 4)`` matrix of candidate critical lambdas for ``segment``.
250 The first ``n`` rows are the box events (a free weight reaching a bound, a
251 blocked multiplier changing sign); the trailing ``p`` rows are the
252 inequality-row events (an inactive row's slack reaching zero, an active
253 row's multiplier changing sign). The generic tracer treats the two blocks
254 uniformly; ``step`` decodes a row index ``>= n`` as a row event.
256 ``state`` is part of the uniform ``ParametricProblem`` signature; the CLA
257 does not need it here because ``segment`` already bundles the active-set
258 masks derived from it.
259 """
260 box = event_ratios(
261 segment.r_alpha,
262 segment.r_beta,
263 segment.gamma,
264 segment.delta,
265 segment.free_in,
266 segment.at_upper,
267 segment.at_lower,
268 self.lower_bounds,
269 self.upper_bounds,
270 )
271 ineq = ineq_event_ratios(
272 segment.r_alpha,
273 segment.r_beta,
274 segment.eta_alpha,
275 segment.eta_beta,
276 segment.active_ineq,
277 self.g_matrix,
278 self.h_vector,
279 )
280 return np.vstack([box, ineq])
282 def step(self, state: TurningPoint, segment: _Segment, sec: int, direction: int, lam: float) -> TurningPoint:
283 """Emit the turning point at ``lam`` after flipping the activity at ``sec``.
285 ``sec < n`` is a box event on asset ``sec``: a "leaves a bound" event
286 (``direction`` in {2, 3}) makes it free, a "moves to a bound" event
287 (``direction`` in {0, 1}) blocks it. ``sec >= n`` is an inequality-row
288 event on row ``sec - n``: ``direction == 0`` activates the row (its slack
289 reached zero), ``direction == 1`` releases it (its multiplier reached
290 zero). The weight vector is continuous across either event.
291 """
292 n = self.dimension
293 free = state.free
294 active_ineq = state.active_ineq
295 if sec < n:
296 free = free.copy()
297 free[sec] = direction >= 2
298 else:
299 active_ineq = active_ineq.copy()
300 active_ineq[sec - n] = direction == 0
301 self._emit(lam, segment.r_alpha + lam * segment.r_beta, free, active_ineq)
302 return self.turning_points[-1]
304 def finish(self, state: TurningPoint, segment: _Segment) -> None:
305 """Emit the minimum-variance endpoint at ``lambda = 0``."""
306 self._emit(0.0, segment.r_alpha, state.free, state.active_ineq)
308 def __len__(self) -> int:
309 """Get the number of turning points in the efficient frontier.
311 Returns:
312 The number of turning points currently stored in the object.
314 """
315 return len(self.turning_points)
317 def _first_turning_point(self) -> TurningPoint:
318 """Calculate the first turning point on the efficient frontier.
320 The first turning point is the maximum-return vertex of the feasible
321 polytope. For the all-ones budget constraint with no inequality rows it is
322 found by the greedy fill of ``init_algo``; for a general equality system
323 ``A w = b`` or any ``G w <= h`` it is found by solving the linear program
324 in ``first_vertex_lp``, which also reports the initially-active rows.
326 Returns:
327 A TurningPoint object representing the first point on the efficient frontier.
329 """
330 if self.g_matrix.shape[0] == 0 and self.a.shape[0] == 1 and np.allclose(self.a, 1.0):
331 return init_algo(
332 mean=self.mean,
333 lower_bounds=self.lower_bounds,
334 upper_bounds=self.upper_bounds,
335 total=float(self.b[0]),
336 )
337 return first_vertex_lp(
338 mean=self.mean,
339 lower_bounds=self.lower_bounds,
340 upper_bounds=self.upper_bounds,
341 a=self.a,
342 b=self.b,
343 tol=self.tol,
344 g=self.g_matrix,
345 h=self.h_vector,
346 )
348 def _append(self, tp: TurningPoint, tol: float | None = None) -> None:
349 """Append a turning point to the list of turning points.
351 This method validates that the turning point satisfies the constraints
352 before adding it to the list.
354 Args:
355 tp: The turning point to append.
356 tol: Tolerance for constraint validation. If None, uses the class's
357 tol attribute. Pass 0 for exact validation.
359 Raises:
360 ValueError: If the turning point violates any constraints.
362 """
363 tol = self.tol if tol is None else tol
365 # (constraint holds?, message if it does not). An empty ``g_matrix`` makes
366 # the inequality ``np.all`` vacuously true, so it never fires when absent.
367 checks: tuple[tuple[bool, str], ...] = (
368 (bool(np.all(tp.weights >= (self.lower_bounds - tol))), "Weights below lower bounds"), # pragma: no mutate
369 (bool(np.all(tp.weights <= (self.upper_bounds + tol))), "Weights above upper bounds"), # pragma: no mutate
370 (
371 bool(np.allclose(self.a @ tp.weights, self.b, atol=1e-7)),
372 "Weights violate the equality constraint A w = b",
373 ),
374 (
375 bool(np.all(self.g_matrix @ tp.weights <= self.h_vector + tol)),
376 "Weights violate the inequality constraint G w <= h",
377 ),
378 )
379 for ok, message in checks:
380 if not ok:
381 raise ValueError(message)
383 self.turning_points.append(tp)
385 def _emit(
386 self,
387 lamb: float,
388 weights: NDArray[np.float64],
389 free: NDArray[np.bool_],
390 active_ineq: NDArray[np.bool_],
391 ) -> None:
392 """Build and store a turning point, projecting away sub-tolerance round-off.
394 Orchestrates the three steps taken at every turning point: refuse the point
395 if the free-asset block is numerically singular (see
396 :meth:`_guard_degeneracy`); project the candidate back onto the feasible
397 set to clear sub-tolerance round-off (see
398 :func:`cvxcla._projection.project_feasible`); then validate and store it
399 (see :meth:`_append`).
401 On tie-heavy or near-degenerate problems (a short, near-rank-deficient
402 sample covariance, duplicated assets, or many coincident events) accumulated
403 floating-point round-off over the many turning points of a large trace can
404 place a free weight a hair outside its box. The covariance there has
405 near-flat directions (its small eigenvalues) and the round-off lies in
406 exactly those directions, so the candidate is optimal to solver precision
407 but not exactly feasible; the projection clears it and is a strict no-op for
408 the well-posed turning points that are already feasible.
409 """
410 self._guard_degeneracy(lamb, free)
411 weights = project_feasible(
412 weights,
413 self.lower_bounds,
414 self.upper_bounds,
415 self.a,
416 self.b,
417 self.g_matrix,
418 self.h_vector,
419 active_ineq,
420 )
421 self._append(TurningPoint(lamb=lamb, weights=weights, free=free, active_ineq=active_ineq))
423 def _guard_degeneracy(self, lamb: float, free: NDArray[np.bool_]) -> None:
424 """Refuse the turning point when the free-asset block is numerically singular.
426 We distinguish two regimes by the conditioning of the free-asset block.
427 While that block stays numerically full rank its solve is reliable and any
428 box violation is round-off, which
429 :func:`cvxcla._projection.project_feasible` clears. Once the
430 free set grows past the covariance rank the block is numerically singular
431 and its solve is unreliable; whatever weights it produces (feasible or not)
432 cannot be trusted, so we refuse and raise an actionable diagnosis instead of
433 silently returning a possibly-suboptimal frontier.
435 The discriminator is the free block's reciprocal condition number, read
436 from its symmetric eigenvalues. Unlike the magnitude of the box violation,
437 which is the residual of a singular solve and therefore varies with the
438 BLAS/LAPACK build, the conditioning is deterministic and portable, so the
439 completed-vs-declined boundary is the same on every platform.
441 The per-turning-point conditioning check is skipped entirely when the full
442 covariance is well conditioned: by interlacing no free block can then be
443 singular, so the check is provably redundant (see
444 :attr:`_free_blocks_well_conditioned`). It runs only when the full
445 covariance is itself near-singular, which is exactly the regime that can
446 produce an untrustworthy free-block solve.
448 Args:
449 lamb: Lambda value of the candidate turning point, used in the message.
450 free: Boolean mask of the free assets at the candidate.
452 Raises:
453 ValueError: With a degeneracy-specific message when the free-asset
454 block is numerically singular (an unreliable solve); otherwise
455 returns without effect.
456 """
457 # When the full covariance clears the floor, interlacing guarantees every
458 # free block does too, so the per-step guard can never fire -- skip it and
459 # the costly per-step rcond. Only a near-singular full covariance needs the
460 # check, and there it runs exactly as before.
461 if not self._free_blocks_well_conditioned:
462 rcond = self.covariance_operator.rcond_free(np.flatnonzero(free))
463 if rcond < _RCOND_FLOOR:
464 n_free = int(np.count_nonzero(free))
465 msg = (
466 f"Critical Line Algorithm hit a degeneracy at lambda={lamb:.4g} "
467 f"(free-set size {n_free}): the free-asset covariance block is "
468 f"numerically singular (reciprocal condition number {rcond:.2g}), "
469 "so its solve is unreliable and the turning point cannot be "
470 "trusted. The trace was stopped rather than risk silently "
471 "returning a suboptimal frontier. This happens when the free set "
472 "grows past the covariance rank (for example a sample covariance "
473 "from far fewer days than assets). Use a well-conditioned, "
474 "full-rank estimate (ample history), or a FactorCovariance backend "
475 "(diagonal-plus-low-rank), which is positive definite by construction."
476 )
477 raise ValueError(msg)
479 @property
480 def frontier(self) -> Frontier:
481 """Get the efficient frontier constructed from the turning points.
483 This property creates a Frontier object from the list of turning points,
484 which can be used to analyze the risk-return characteristics of the
485 efficient portfolios.
487 Returns:
488 A Frontier object representing the efficient frontier.
490 """
491 return Frontier(
492 covariance=self.covariance,
493 mean=self.mean,
494 frontier=[FrontierPoint(point.weights) for point in self.turning_points],
495 )
498class ProblemBuilder:
499 """Chainable builder that assembles the polyhedral pieces of a CLA problem.
501 A thin, chainable convenience layer over the explicit :class:`CLA`
502 constructor. It exists purely for readability: portfolio practitioners expect
503 to say "long-only, fully invested" rather than to remember that the budget is
504 encoded as ``a=np.ones((1, n)), b=np.ones(1)``. Every method maps one-to-one
505 onto a constructor argument, so the builder adds no modelling power and imposes
506 no expression algebra: it accepts the same polyhedral pieces the CLA already
507 supports (a quadratic objective, box bounds, linear equalities ``A w = b``, and
508 linear inequalities ``G w <= h``) and nothing else. Anything the explicit
509 constructor cannot trace, the builder cannot express either.
511 Construct one via :meth:`CLA.problem`, chain the constraint methods (each
512 returns ``self``), and finish with :meth:`trace`, which builds the ``CLA`` and
513 runs the full parametric trace, returning the solved object whose ``frontier``
514 and ``turning_points`` describe the entire efficient frontier (not a single
515 optimum, which is the distinction from a one-shot convex solver).
517 Attributes:
518 mean: Vector of expected returns, fixing the problem dimension ``n``.
519 covariance: The covariance, either a plain ``numpy`` array or a
520 ``QuadraticForm`` backend (e.g. ``FactorCovariance``), passed through
521 to ``CLA`` unchanged so the structured backends keep their advantage.
523 Examples:
524 >>> import numpy as np
525 >>> from cvxcla import CLA
526 >>> rng = np.random.default_rng(0)
527 >>> mean = rng.uniform(0.0, 1.0, 4)
528 >>> covariance = np.eye(4)
529 >>> cla = CLA.problem(mean, covariance).long_only().budget().trace()
530 >>> len(cla) > 0
531 True
532 """
534 def __init__(self, mean: NDArray[np.float64], covariance: NDArray[np.float64] | QuadraticForm) -> None:
535 """Start a builder for an ``n``-asset problem.
537 Args:
538 mean: Vector of expected returns of length ``n``.
539 covariance: Covariance matrix or ``QuadraticForm`` backend.
540 """
541 self.mean = np.asarray(mean, dtype=np.float64)
542 self.covariance = covariance
543 self._lower: NDArray[np.float64] | None = None
544 self._upper: NDArray[np.float64] | None = None
545 self._a_blocks: list[NDArray[np.float64]] = []
546 self._b_blocks: list[NDArray[np.float64]] = []
547 self._g_blocks: list[NDArray[np.float64]] = []
548 self._h_blocks: list[NDArray[np.float64]] = []
550 @property
551 def _n(self) -> int:
552 """Number of assets ``n``, fixed by ``mean``."""
553 return int(self.mean.shape[0])
555 def _as_vector(self, value: float | NDArray[np.float64], name: str) -> NDArray[np.float64]:
556 """Broadcast a scalar or length-``n`` array to a length-``n`` vector.
558 Args:
559 value: A scalar (applied to every asset) or a length-``n`` array.
560 name: Argument name, used in the error message.
562 Returns:
563 A fresh length-``n`` float array.
565 Raises:
566 ValueError: If an array is passed whose length is not ``n``.
567 """
568 array = np.asarray(value, dtype=np.float64)
569 if array.ndim == 0:
570 return np.full(self._n, float(array))
571 if array.shape != (self._n,):
572 msg = f"{name} must be a scalar or a length-{self._n} vector, got shape {array.shape}"
573 raise ValueError(msg)
574 return array.astype(np.float64, copy=True)
576 def bounds(self, lower: float | NDArray[np.float64], upper: float | NDArray[np.float64]) -> "ProblemBuilder":
577 """Set the box bounds ``lower <= w <= upper``.
579 Args:
580 lower: Lower bound, a scalar (same for every asset) or length-``n`` array.
581 upper: Upper bound, a scalar or length-``n`` array.
583 Returns:
584 ``self``, for chaining.
585 """
586 self._lower = self._as_vector(lower, "lower")
587 self._upper = self._as_vector(upper, "upper")
588 return self
590 def long_only(self, upper: float | NDArray[np.float64] = 1.0) -> "ProblemBuilder":
591 """Set long-only box bounds ``0 <= w <= upper`` (``upper`` defaults to ``1``).
593 Args:
594 upper: Upper bound, a scalar or length-``n`` array; defaults to ``1.0``.
596 Returns:
597 ``self``, for chaining.
598 """
599 return self.bounds(0.0, upper)
601 def budget(self, total: float = 1.0) -> "ProblemBuilder":
602 """Add the fully-invested budget constraint ``sum(w) = total``.
604 This is the canonical all-ones equality row; ``total=0`` gives a
605 dollar-neutral book. Equivalent to ``equality(np.ones(n), total)``.
607 Args:
608 total: The right-hand side of ``sum(w) = total``; defaults to ``1.0``.
610 Returns:
611 ``self``, for chaining.
612 """
613 return self.equality(np.ones(self._n), total)
615 def equality(self, a: NDArray[np.float64], b: float | NDArray[np.float64]) -> "ProblemBuilder":
616 """Add one or more equality rows ``A w = b``.
618 Accepts a single row (a length-``n`` vector with a scalar right-hand side)
619 or a block of rows (an ``(m, n)`` matrix with a length-``m`` right-hand
620 side). Repeated calls accumulate rows, so a budget plus a sector-neutrality
621 block can be added separately.
623 Args:
624 a: A length-``n`` row vector or an ``(m, n)`` matrix.
625 b: The matching right-hand side: a scalar for a single row, or a
626 length-``m`` vector for a block.
628 Returns:
629 ``self``, for chaining.
631 Raises:
632 ValueError: If ``a`` does not have ``n`` columns, or ``b``'s length
633 does not match the number of rows of ``a``.
634 """
635 a_block = np.atleast_2d(np.asarray(a, dtype=np.float64))
636 b_block = np.atleast_1d(np.asarray(b, dtype=np.float64))
637 self._validate_rows(a_block, b_block, "equality", "b")
638 self._a_blocks.append(a_block)
639 self._b_blocks.append(b_block)
640 return self
642 def inequality(self, g: NDArray[np.float64], h: float | NDArray[np.float64]) -> "ProblemBuilder":
643 """Add one or more inequality rows ``G w <= h``.
645 Like :meth:`equality` but for ``<=`` rows (e.g. a group- or
646 sector-exposure cap). A ``>=`` row is expressed by negating both ``g`` and
647 ``h``. Repeated calls accumulate rows.
649 Args:
650 g: A length-``n`` row vector or a ``(p, n)`` matrix.
651 h: The matching right-hand side: a scalar for a single row, or a
652 length-``p`` vector for a block.
654 Returns:
655 ``self``, for chaining.
657 Raises:
658 ValueError: If ``g`` does not have ``n`` columns, or ``h``'s length
659 does not match the number of rows of ``g``.
660 """
661 g_block = np.atleast_2d(np.asarray(g, dtype=np.float64))
662 h_block = np.atleast_1d(np.asarray(h, dtype=np.float64))
663 self._validate_rows(g_block, h_block, "inequality", "h")
664 self._g_blocks.append(g_block)
665 self._h_blocks.append(h_block)
666 return self
668 def _validate_rows(self, lhs: NDArray[np.float64], rhs: NDArray[np.float64], method: str, rhs_name: str) -> None:
669 """Check a constraint block has ``n`` columns and a matching right-hand side.
671 Args:
672 lhs: The ``(m, n)`` coefficient block.
673 rhs: The length-``m`` right-hand side.
674 method: The calling method name, used in error messages.
675 rhs_name: The right-hand-side argument name, used in error messages.
677 Raises:
678 ValueError: If the column count is not ``n`` or the lengths disagree.
679 """
680 if lhs.shape[1] != self._n:
681 msg = f"{method}: coefficient matrix must have {self._n} columns, got shape {lhs.shape}"
682 raise ValueError(msg)
683 if rhs.shape[0] != lhs.shape[0]:
684 msg = f"{method}: {rhs_name} must have {lhs.shape[0]} entries to match the rows, got {rhs.shape[0]}"
685 raise ValueError(msg)
687 def trace(self) -> CLA:
688 """Assemble the pieces, build the ``CLA``, and run the full trace.
690 Returns:
691 The solved :class:`CLA`, whose ``frontier`` and ``turning_points``
692 describe the entire efficient frontier.
694 Raises:
695 ValueError: If no box bounds were set (call :meth:`bounds` or
696 :meth:`long_only`), or no equality constraint was added (call
697 :meth:`budget` or :meth:`equality`).
698 """
699 lower, upper = self._resolved_bounds()
700 if not self._a_blocks:
701 msg = "a CLA problem needs an equality constraint: call .budget() or .equality(A, b)"
702 raise ValueError(msg)
704 g = np.vstack(self._g_blocks) if self._g_blocks else None
705 h = np.concatenate(self._h_blocks) if self._h_blocks else None
706 return CLA(
707 mean=self.mean,
708 covariance=self.covariance,
709 lower_bounds=lower,
710 upper_bounds=upper,
711 a=np.vstack(self._a_blocks),
712 b=np.concatenate(self._b_blocks),
713 g=g,
714 h=h,
715 )
717 def _resolved_bounds(self) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
718 """Return the box bounds, raising if they were never set.
720 Returns:
721 The ``(lower, upper)`` box-bound vectors.
723 Raises:
724 ValueError: If no box bounds were set (call :meth:`bounds` or
725 :meth:`long_only`).
726 """
727 if self._lower is None or self._upper is None:
728 msg = "set box bounds before tracing: call .long_only() or .bounds(lower, upper)"
729 raise ValueError(msg)
730 return self._lower, self._upper