Coverage for src/cvxcla/cla.py: 100%
123 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-09-29 05:24 +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 cast
14import numpy as np
15from numpy.typing import NDArray
17from ._builders import ProblemBuilder
18from ._checks import check_feasible, guard_degeneracy, well_conditioned
19from ._events import segment_events
20from ._kkt import Segment, critical_segment
21from ._leverage import LeverageLift, SignedLift, mask_leg_events, tighten_at_minimum_gross
22from ._projection import project_feasible
23from .first import classify_vertex, first_turning_point
24from .operators import DenseCovariance, QuadraticForm
25from .pathtracer import InequalityConstrained, trace
26from .types import Frontier, FrontierPoint, TurningPoint
29@dataclass(frozen=True)
30class CLA(InequalityConstrained):
31 """Critical Line Algorithm implementation based on Markowitz's approach.
33 This class implements the Critical Line Algorithm as described by Harry Markowitz
34 and colleagues. It computes the entire efficient frontier by finding all turning
35 points, which are the points where the set of assets at their bounds changes.
37 The algorithm starts with the first turning point (the portfolio with the highest
38 expected return) and then iteratively computes the next turning point with a lower
39 expected return until it reaches the minimum variance portfolio.
41 Attributes:
42 mean: Vector of expected returns for each asset.
43 covariance: Covariance matrix of asset returns, either as a plain
44 ``numpy`` array or as a ``CovarianceOperator`` backend
45 (see ``cvxcla.operators``).
46 lower_bounds: Vector of lower bounds for asset weights.
47 upper_bounds: Vector of upper bounds for asset weights.
48 a: Equality-constraint matrix ``A`` of ``A w = b`` (``m x n``). The
49 canonical case is the single all-ones budget row (``sum(w) = b``),
50 but an arbitrary equality system is supported: weighted single rows
51 and ``m > 1`` rows (e.g. budget plus sector- or factor-neutrality).
52 The all-ones budget (any right-hand side, including ``0`` for
53 dollar-neutral) uses the greedy first vertex of
54 :func:`cvxcla.first.init_algo`; a general ``A`` uses the
55 linear-programming first vertex of
56 :func:`cvxcla.first.first_vertex_lp`.
57 b: Equality-constraint right-hand side ``b`` (length ``m``); ``[1]`` for
58 the fully-invested budget, ``[0]`` for dollar-neutral, and so on.
59 g: Optional inequality-constraint matrix ``G`` of ``G w <= h``
60 (``p x n``), e.g. a group- or sector-exposure cap. ``None`` (the
61 default) means no inequality rows, recovering the equality-only
62 problem exactly. A ``>=`` constraint is expressed by negating both
63 ``g`` and ``h``. Each *active* row (held at equality ``g_i w = h_i``)
64 enters the reduced KKT system as an extra equality row, so the
65 covariance is still touched only through the ``QuadraticForm``
66 interface; box bounds remain a separate per-variable active set.
67 h: Optional inequality-constraint right-hand side ``h`` (length ``p``).
68 turning_points: List of turning points on the efficient frontier.
69 tol: Tolerance for numerical calculations.
70 logger: Logger instance for logging information and errors.
71 leverage: Optional cap ``c`` on the gross exposure, ``||w||_1 <= c``.
72 ``None`` (the default) means no cap. The 1-norm is traced exactly by
73 splitting every asset whose box straddles zero into a long and a short
74 leg (see :mod:`cvxcla._leverage`); the turning points are reported in
75 the original asset weights, and ``active_ineq`` covers the rows of
76 ``g`` only.
78 """
80 mean: NDArray[np.float64]
81 covariance: NDArray[np.float64] | QuadraticForm
82 lower_bounds: NDArray[np.float64]
83 upper_bounds: NDArray[np.float64]
84 a: NDArray[np.float64]
85 b: NDArray[np.float64]
86 g: NDArray[np.float64] | None = None
87 h: NDArray[np.float64] | None = None
88 turning_points: list[TurningPoint] = field(default_factory=list)
89 tol: float = 1e-5 # pragma: no mutate
90 logger: logging.Logger = field(default_factory=lambda: logging.getLogger(__name__))
91 leverage: float | None = None
93 @classmethod
94 def problem(
95 cls, mean: NDArray[np.float64], covariance: NDArray[np.float64] | QuadraticForm
96 ) -> "ProblemBuilder[CLA]":
97 """Start a fluent :class:`cvxcla._builders.ProblemBuilder` for this problem.
99 A readability convenience over the explicit constructor: chain
100 ``.long_only()``/``.budget()``/``.equality()``/``.inequality()`` and finish
101 with ``.trace()``. The builder maps one-to-one onto the constructor
102 arguments and adds no modelling power.
104 ``cls`` is handed to the builder as the class it should construct, which is
105 what lets the builder live in a leaf module that never imports this one.
107 Args:
108 mean: Vector of expected returns of length ``n``.
109 covariance: Covariance matrix or ``QuadraticForm`` backend.
111 Returns:
112 A ``ProblemBuilder`` ready to accept constraints, whose ``.trace()``
113 returns a solved ``CLA``.
114 """
115 return ProblemBuilder(mean, covariance, solver=cls)
117 @cached_property
118 def covariance_operator(self) -> QuadraticForm:
119 """Return the covariance as a ``QuadraticForm`` backend.
121 A plain ``numpy`` covariance matrix is wrapped in ``DenseCovariance``;
122 an object already implementing the protocol is passed through. This is
123 the single point where the input form is normalised.
124 """
125 if isinstance(self.covariance, QuadraticForm):
126 return self.covariance
127 return DenseCovariance(self.covariance)
129 @property
130 def dimension(self) -> int:
131 """Number of assets ``n`` (the problem dimension for the path tracer)."""
132 return len(self.mean)
134 @cached_property
135 def _free_blocks_well_conditioned(self) -> bool:
136 """Whether every free-block solve is numerically safe (see :func:`cvxcla._checks.well_conditioned`).
138 Decided once, up front: when ``True`` the per-turning-point guard in
139 :meth:`_emit` is provably redundant and skipped.
140 """
141 return well_conditioned(self.covariance_operator)
143 def __post_init__(self) -> None:
144 """Initialize the CLA object and compute the efficient frontier.
146 This method is automatically called after initialization. It computes
147 the entire efficient frontier by finding all turning points, starting
148 from the first turning point (highest expected return) and iteratively
149 computing the next turning point with a lower expected return until
150 it reaches the minimum variance portfolio.
152 The actual walk is driven by the generic ``cvxcla.pathtracer.trace``
153 loop; this class supplies the portfolio-specific hooks (``begin``,
154 ``segment``, ``event_matrix``, ``step``, ``finish``) it calls.
156 The reduced KKT system at each turning point is solved by block
157 elimination: a single multi-RHS solve against the free covariance block
158 (via the covariance backend), covering the constraint columns and the
159 alpha and beta systems together so ``Sigma_FF`` is factorised once, then a
160 small Schur-complement solve ``A_F @ Sigma_FF^{-1} @ A_F.T`` over the
161 equality (and active inequality) rows. The covariance only enters through
162 the ``QuadraticForm`` interface, so structured backends (e.g.
163 ``FactorCovariance``) never materialise an n x n matrix.
165 Raises:
166 RuntimeError: If all variables are blocked, which would make the
167 system of equations singular.
168 ValueError: If the inequality matrix ``g`` and vector ``h`` have
169 mismatched or wrong shapes, or ``leverage`` is not a
170 positive finite number.
172 """
173 if self.g_matrix.shape[1] != self.dimension:
174 msg = f"g must have {self.dimension} columns, got shape {self.g_matrix.shape}"
175 raise ValueError(msg)
176 if self.h_vector.shape[0] != self.g_matrix.shape[0]:
177 msg = f"h must have {self.g_matrix.shape[0]} entries, got {self.h_vector.shape[0]}"
178 raise ValueError(msg)
179 if self.leverage is None:
180 trace(self)
181 return
182 if not (np.isfinite(self.leverage) and self.leverage > 0):
183 msg = f"leverage must be a positive finite number, got {self.leverage}"
184 raise ValueError(msg)
185 self._trace_leveraged(self.leverage)
187 def _trace_leveraged(self, leverage: float) -> None:
188 """Trace the frontier under ``||w||_1 <= leverage`` via the signed lift.
190 Builds the lifted problem over the long/short legs -- the covariance as a
191 :class:`cvxcla._leverage.SignedLift` of this problem's backend, every
192 constraint matrix mapped through ``w = P x``, and the gross-exposure row
193 ``sum(x) <= leverage`` appended to ``G`` -- traces it with
194 :class:`_LeveragedCLA`, and maps its turning points back to asset weights.
195 A cap at the smallest feasible gross exposure is first resolved into
196 tightened bounds (see :func:`cvxcla._leverage.tighten_at_minimum_gross`),
197 which keeps the maximum-return vertex non-degenerate.
198 The objective ``x.T P.T Sigma P x / 2 - lam mean.T P x`` is the original one
199 in ``w = P x``, so the lifted ``lambda`` is the original ``lambda``.
201 Args:
202 leverage: The gross-exposure cap ``c``.
203 """
204 lower, upper, keep_cap = tighten_at_minimum_gross(
205 self.lower_bounds, self.upper_bounds, self.a, self.b, self.g_matrix, self.h_vector, leverage, self.tol
206 )
207 lift = LeverageLift.from_bounds(lower, upper)
208 g, h = lift.with_cap(self.g_matrix, self.h_vector, leverage if keep_cap else None)
209 lifted = _LeveragedCLA(
210 mean=self.mean[lift.asset] * lift.sign,
211 covariance=SignedLift(self.covariance_operator, lift.asset, lift.sign),
212 lower_bounds=lift.lower,
213 upper_bounds=lift.upper,
214 a=lift.columns(self.a),
215 b=self.b,
216 g=g,
217 h=h,
218 tol=self.tol,
219 logger=self.logger,
220 lift=lift,
221 )
222 p = self.g_matrix.shape[0]
223 for tp in lifted.turning_points:
224 self._append(lift.to_turning_point(tp, self.dimension, p))
226 def begin(self) -> tuple[float, TurningPoint]:
227 """Record the first turning point and start the trace at ``lambda = inf``.
229 Returns:
230 ``(inf, first_turning_point)``: the starting lambda bound and the
231 initial state for the path tracer.
232 """
233 first = self._first_turning_point()
234 self._append(first)
235 return np.inf, first
237 def segment(self, state: TurningPoint) -> Segment:
238 """Solve the reduced KKT system for the critical-line segment at ``state``."""
239 return critical_segment(
240 self.covariance_operator,
241 self.mean,
242 self.a,
243 self.b,
244 self.g_matrix,
245 self.h_vector,
246 self.lower_bounds,
247 self.upper_bounds,
248 self.tol,
249 state,
250 )
252 def event_matrix(self, state: TurningPoint, segment: Segment) -> NDArray[np.float64]: # noqa: ARG002
253 """Return the ``(n + p, 4)`` event matrix for ``segment`` (see :func:`cvxcla._events.segment_events`).
255 The first ``n`` rows are box events and the trailing ``p`` rows are
256 inequality-row events; ``step`` decodes a row index ``>= n`` as a row
257 event. ``state`` is part of the uniform ``ParametricProblem`` signature;
258 the CLA does not need it because ``segment`` already bundles the masks.
259 """
260 return segment_events(segment, self.lower_bounds, self.upper_bounds, self.g_matrix, self.h_vector)
262 def step(self, state: TurningPoint, segment: Segment, sec: int, direction: int, lam: float) -> TurningPoint:
263 """Emit the turning point at ``lam`` after flipping the activity at ``sec``.
265 ``sec < n`` is a box event on asset ``sec``: a "leaves a bound" event
266 (``direction`` in {2, 3}) makes it free, a "moves to a bound" event
267 (``direction`` in {0, 1}) blocks it. ``sec >= n`` is an inequality-row
268 event on row ``sec - n``: ``direction == 0`` activates the row (its slack
269 reached zero), ``direction == 1`` releases it (its multiplier reached
270 zero). The weight vector is continuous across either event.
271 """
272 n = self.dimension
273 free = state.free
274 active_ineq = state.active_ineq
275 if sec < n:
276 free = free.copy()
277 free[sec] = direction >= 2
278 else:
279 active_ineq = active_ineq.copy()
280 active_ineq[sec - n] = direction == 0
281 self._emit(lam, segment.r_alpha + lam * segment.r_beta, free, active_ineq)
282 return self.turning_points[-1]
284 def finish(self, state: TurningPoint, segment: Segment) -> None:
285 """Emit the minimum-variance endpoint at ``lambda = 0``."""
286 self._emit(0.0, segment.r_alpha, state.free, state.active_ineq)
288 def __len__(self) -> int:
289 """Get the number of turning points in the efficient frontier.
291 Returns:
292 The number of turning points currently stored in the object.
294 """
295 return len(self.turning_points)
297 def _first_turning_point(self) -> TurningPoint:
298 """Return the maximum-return vertex (see :func:`cvxcla.first.first_turning_point`)."""
299 return first_turning_point(
300 self.mean, self.lower_bounds, self.upper_bounds, self.a, self.b, self.g_matrix, self.h_vector, self.tol
301 )
303 def _append(self, tp: TurningPoint, tol: float | None = None) -> None:
304 """Append a turning point to the list of turning points.
306 This method validates that the turning point satisfies the constraints
307 (see :func:`cvxcla._checks.check_feasible`) before adding it to the list.
309 Args:
310 tp: The turning point to append.
311 tol: Tolerance for constraint validation. If None, uses the class's
312 tol attribute. Pass 0 for exact validation.
314 Raises:
315 ValueError: If the turning point violates any constraints.
317 """
318 tol = self.tol if tol is None else tol
319 check_feasible(
320 tp.weights,
321 self.lower_bounds,
322 self.upper_bounds,
323 self.a,
324 self.b,
325 self.g_matrix,
326 self.h_vector,
327 self.leverage,
328 tol,
329 )
330 self.turning_points.append(tp)
332 def _emit(
333 self,
334 lamb: float,
335 weights: NDArray[np.float64],
336 free: NDArray[np.bool_],
337 active_ineq: NDArray[np.bool_],
338 ) -> None:
339 """Build and store a turning point, projecting away sub-tolerance round-off.
341 Orchestrates the three steps taken at every turning point: refuse the point
342 if the free-asset block is numerically singular (see
343 :func:`cvxcla._checks.guard_degeneracy`); project the candidate back onto the feasible
344 set to clear sub-tolerance round-off (see
345 :func:`cvxcla._projection.project_feasible`); then validate and store it
346 (see :meth:`_append`).
348 On tie-heavy or near-degenerate problems (a short, near-rank-deficient
349 sample covariance, duplicated assets, or many coincident events) accumulated
350 floating-point round-off over the many turning points of a large trace can
351 place a free weight a hair outside its box. The covariance there has
352 near-flat directions (its small eigenvalues) and the round-off lies in
353 exactly those directions, so the candidate is optimal to solver precision
354 but not exactly feasible; the projection clears it and is a strict no-op for
355 the well-posed turning points that are already feasible.
356 """
357 if not self._free_blocks_well_conditioned:
358 guard_degeneracy(self.covariance_operator, lamb, free)
359 weights = project_feasible(
360 weights,
361 self.lower_bounds,
362 self.upper_bounds,
363 self.a,
364 self.b,
365 self.g_matrix,
366 self.h_vector,
367 active_ineq,
368 )
369 self._append(TurningPoint(lamb=lamb, weights=weights, free=free, active_ineq=active_ineq))
371 @property
372 def frontier(self) -> Frontier:
373 """Get the efficient frontier constructed from the turning points.
375 This property creates a Frontier object from the list of turning points,
376 which can be used to analyze the risk-return characteristics of the
377 efficient portfolios.
379 Returns:
380 A Frontier object representing the efficient frontier.
382 """
383 return Frontier(
384 covariance=self.covariance,
385 mean=self.mean,
386 frontier=[FrontierPoint(point.weights) for point in self.turning_points],
387 )
390@dataclass(frozen=True)
391class _LeveragedCLA(CLA):
392 """The CLA over the long/short legs of a leverage-constrained problem.
394 Built only by :meth:`CLA._trace_leveraged`; its covariance is a
395 :class:`cvxcla._leverage.SignedLift` and its last inequality row is the
396 gross-exposure cap. It differs from the plain CLA in three places, all
397 stemming from the lifted covariance being singular along "raise both legs of
398 one asset":
400 * the leave-a-bound event of a leg whose partner is off its lower bound is
401 masked, so both legs of an asset are never free together (see
402 :func:`cvxcla._leverage.mask_leg_events`);
403 * a maximum-return vertex with overlapping legs is netted before the trace
404 starts (the LP can return one when the cap is tight with a zero multiplier);
405 * the up-front conditioning test is taken on the asset covariance, since the
406 lifted form is singular as a whole but never on a free block that is traced.
408 Attributes:
409 lift: The signed leg structure.
410 """
412 lift: LeverageLift | None = None
414 @property
415 def _legs(self) -> LeverageLift:
416 """The leg structure (always set by :meth:`CLA._trace_leveraged`)."""
417 return cast(LeverageLift, self.lift)
419 @cached_property
420 def _free_blocks_well_conditioned(self) -> bool:
421 """Whether the asset covariance clears the singularity floor.
423 A traced free block holds at most one leg per asset, so it is a signed
424 principal block of the asset covariance and interlacing applies to that.
425 """
426 return well_conditioned(cast(SignedLift, self.covariance).base)
428 def event_matrix(self, state: TurningPoint, segment: Segment) -> NDArray[np.float64]:
429 """Return the event matrix with the competing leg events masked."""
430 events = super().event_matrix(state, segment)
431 legs = self.dimension
432 events[:legs] = mask_leg_events(events[:legs], self._legs.partner, segment.at_lower)
433 return events
435 def _first_turning_point(self) -> TurningPoint:
436 """Return the maximum-return vertex with overlapping legs netted out."""
437 first = super()._first_turning_point()
438 netted = self._legs.net(first.weights)
439 if np.array_equal(netted, first.weights):
440 return first
441 return classify_vertex(
442 netted, self.lower_bounds, self.upper_bounds, self.a, self.g_matrix, self.h_vector, self.tol
443 )