Coverage for src/cvxcla/first.py: 100%
54 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-09-29 05:24 +0000
1"""First turning point computation for the Critical Line Algorithm.
3This module provides functions to compute the first turning point on the efficient frontier,
4which is the portfolio with the highest expected return that satisfies the constraints.
5Two implementations are provided: a direct algorithm and a linear programming approach.
6"""
8from __future__ import annotations
10import numpy as np
11from numpy.typing import NDArray
12from scipy.optimize import linprog # type: ignore[import-untyped]
14from .types import TurningPoint
17#
18def init_algo(
19 mean: NDArray[np.float64],
20 lower_bounds: NDArray[np.float64],
21 upper_bounds: NDArray[np.float64],
22 total: float = 1.0,
23) -> TurningPoint:
24 """Compute the first turning point for a single all-ones budget constraint.
26 The key insight behind Markowitz's CLA is to find first the
27 turning point associated with the highest expected return, and then
28 compute the sequence of turning points, each with a lower expected
29 return than the previous. That first turning point consists in the
30 smallest subset of assets with highest return such that the sum of
31 their upper boundaries equals or exceeds the budget ``total``.
33 We sort the expected returns in descending order.
34 This gives us a sequence for searching for the
35 first free asset. All weights are initially set to their lower bounds,
36 and following the sequence from the previous step, we move those
37 weights from the lower to the upper bound until the sum of weights
38 reaches ``total``. The last iterated weight is then reduced
39 to comply with the constraint that the sum of weights equals ``total``.
40 This last weight is the first free asset,
41 and the resulting vector of weights the first turning point.
43 Args:
44 mean: Vector of expected returns.
45 lower_bounds: Lower box bounds.
46 upper_bounds: Upper box bounds.
47 total: Target sum of weights (the right-hand side ``b`` of the all-ones
48 budget constraint ``sum(w) = total``; ``1`` for fully-invested,
49 ``0`` for dollar-neutral, ``> 1`` for a leveraged total).
50 """
51 if np.any(lower_bounds > upper_bounds):
52 msg = "Lower bounds must be less than or equal to upper bounds"
53 raise ValueError(msg)
55 # Initialize weights to lower bounds
56 weights = np.copy(lower_bounds).astype(np.float64)
57 free = np.full_like(mean, False, dtype=np.bool_)
59 # Move weights from lower to upper bound until the sum reaches ``total``. The
60 # check needs a tolerance: the increment ``total - sum(weights)`` can bring the
61 # sum to ``total`` only up to floating-point error, and without the slack the
62 # loop would move on and mark the NEXT asset (sitting on its bound) as free
63 # while the genuinely interior asset stays blocked.
64 for index in np.argsort(-mean):
65 weights[index] += np.min([upper_bounds[index] - lower_bounds[index], total - np.sum(weights)])
66 if np.sum(weights) >= total - 1e-12:
67 free[index] = True
68 break
70 if not np.any(free):
71 # No asset ended up interior: the bounds cannot sum to the target.
72 msg = "Could not construct a fully invested portfolio"
73 raise ValueError(msg)
75 # Return first turning point, the point with the highest expected return.
76 return TurningPoint(free=free, weights=weights)
79def first_vertex_lp(
80 mean: NDArray[np.float64],
81 lower_bounds: NDArray[np.float64],
82 upper_bounds: NDArray[np.float64],
83 a: NDArray[np.float64],
84 b: NDArray[np.float64],
85 tol: float,
86 g: NDArray[np.float64] | None = None,
87 h: NDArray[np.float64] | None = None,
88) -> TurningPoint:
89 """Compute the first turning point for a general ``A w = b``, ``G w <= h`` system.
91 The maximum-return vertex of the feasible polytope
92 ``{w : A w = b, G w <= h, lower <= w <= upper}`` is a linear program,
93 ``maximize mean @ w``. The greedy fill of :func:`init_algo` only solves the
94 single all-ones budget with no inequality rows; for a general (weighted, or
95 multi-row) ``A`` or any ``G`` we solve the LP directly with HiGHS (via
96 :func:`scipy.optimize.linprog`), which returns a vertex. The free set is read
97 off the solution (assets strictly inside their box bounds) and the initial
98 active inequality set off the tight rows (``g_i w`` at ``h_i`` to tolerance).
100 Args:
101 mean: Vector of expected returns.
102 lower_bounds: Lower box bounds.
103 upper_bounds: Upper box bounds.
104 a: Equality-constraint matrix (``m x n``).
105 b: Equality-constraint right-hand side (length ``m``).
106 tol: Tolerance for classifying an asset as free (strictly interior) and a
107 row as active (tight).
108 g: Inequality-constraint matrix (``p x n``); ``None`` means no rows.
109 h: Inequality-constraint right-hand side (length ``p``).
111 Returns:
112 The maximum-return vertex as a :class:`TurningPoint`, carrying the active
113 inequality rows in ``active_ineq``.
115 Raises:
116 ValueError: If the linear program is infeasible or unbounded (the
117 constraints admit no maximum-return vertex), or if that vertex is
118 degenerate: the free set does not span the equality rows together with
119 the active inequality rows, so the reduced KKT system would be
120 singular. That case is declined here rather than left to surface as an
121 opaque singular-matrix error later in the trace.
122 """
123 g = np.zeros((0, mean.shape[0])) if g is None else np.asarray(g, dtype=np.float64)
124 h = np.zeros(0) if h is None else np.asarray(h, dtype=np.float64)
126 weights = _solve_max_return_lp(mean, lower_bounds, upper_bounds, a, b, g, h)
127 return classify_vertex(weights, lower_bounds, upper_bounds, a, g, h, tol)
130def first_turning_point(
131 mean: NDArray[np.float64],
132 lower_bounds: NDArray[np.float64],
133 upper_bounds: NDArray[np.float64],
134 a: NDArray[np.float64],
135 b: NDArray[np.float64],
136 g: NDArray[np.float64],
137 h: NDArray[np.float64],
138 tol: float,
139) -> TurningPoint:
140 """Calculate the first turning point on the efficient frontier.
142 The first turning point is the maximum-return vertex of the feasible
143 polytope. For the all-ones budget constraint with no inequality rows it is
144 found by the greedy fill of :func:`init_algo`; for a general equality system
145 ``A w = b`` or any ``G w <= h`` it is found by solving the linear program
146 in :func:`first_vertex_lp`, which also reports the initially-active rows.
148 Args:
149 mean: Vector of expected returns.
150 lower_bounds: Lower box bounds.
151 upper_bounds: Upper box bounds.
152 a: Equality-constraint matrix ``A`` of ``A w = b``.
153 b: Equality-constraint right-hand side ``b``.
154 g: Inequality-constraint matrix ``G`` of ``G w <= h`` (``(p, n)``).
155 h: Inequality-constraint right-hand side ``h`` (length ``p``).
156 tol: Tolerance for the linear-programming vertex classification.
158 Returns:
159 A TurningPoint object representing the first point on the efficient frontier.
160 """
161 if g.shape[0] == 0 and a.shape[0] == 1 and np.allclose(a, 1.0):
162 return init_algo(mean=mean, lower_bounds=lower_bounds, upper_bounds=upper_bounds, total=float(b[0]))
163 return first_vertex_lp(mean=mean, lower_bounds=lower_bounds, upper_bounds=upper_bounds, a=a, b=b, tol=tol, g=g, h=h)
166def classify_vertex(
167 weights: NDArray[np.float64],
168 lower_bounds: NDArray[np.float64],
169 upper_bounds: NDArray[np.float64],
170 a: NDArray[np.float64],
171 g: NDArray[np.float64],
172 h: NDArray[np.float64],
173 tol: float,
174) -> TurningPoint:
175 """Read the free set and the active rows off a maximum-return vertex.
177 An asset is free when it sits strictly inside its box (by more than ``tol``)
178 and an inequality row is active when it is tight to ``tol``. The vertex is
179 then checked for degeneracy (see :func:`_reject_degenerate_vertex`).
181 Args:
182 weights: The vertex weights.
183 lower_bounds: Lower box bounds.
184 upper_bounds: Upper box bounds.
185 a: Equality-constraint matrix (``m x n``).
186 g: Inequality-constraint matrix (``p x n``); empty ``(0, n)`` when none.
187 h: Inequality-constraint right-hand side (length ``p``).
188 tol: Classification tolerance.
190 Returns:
191 The vertex as a :class:`TurningPoint` carrying its active rows.
193 Raises:
194 ValueError: If the vertex is degenerate.
195 """
196 free = (weights > lower_bounds + tol) & (weights < upper_bounds - tol)
197 active_ineq = (g @ weights >= h - tol) if g.shape[0] else np.zeros(0, dtype=bool)
199 _reject_degenerate_vertex(a, g, free, active_ineq)
200 return TurningPoint(free=free, weights=weights, active_ineq=active_ineq)
203def _solve_max_return_lp(
204 mean: NDArray[np.float64],
205 lower_bounds: NDArray[np.float64],
206 upper_bounds: NDArray[np.float64],
207 a: NDArray[np.float64],
208 b: NDArray[np.float64],
209 g: NDArray[np.float64],
210 h: NDArray[np.float64],
211) -> NDArray[np.float64]:
212 """Solve the maximum-return linear program and return its vertex weights.
214 ``maximize mean @ w`` (as ``minimize -mean @ w``) subject to ``A w = b``,
215 ``G w <= h`` and the box bounds, via HiGHS. The inequality rows are passed
216 only when ``g`` is non-empty.
218 Args:
219 mean: Vector of expected returns.
220 lower_bounds: Lower box bounds.
221 upper_bounds: Upper box bounds.
222 a: Equality-constraint matrix (``m x n``).
223 b: Equality-constraint right-hand side (length ``m``).
224 g: Inequality-constraint matrix (``p x n``); empty ``(0, n)`` when none.
225 h: Inequality-constraint right-hand side (length ``p``).
227 Returns:
228 The vertex weights ``w`` as a 1d ``float64`` array.
230 Raises:
231 ValueError: If the linear program is infeasible or unbounded.
232 """
233 has_ineq = g.shape[0] > 0
234 result = linprog(
235 c=-np.asarray(mean, dtype=np.float64),
236 A_eq=np.asarray(a, dtype=np.float64),
237 b_eq=np.asarray(b, dtype=np.float64),
238 A_ub=g if has_ineq else None,
239 b_ub=h if has_ineq else None,
240 bounds=list(zip(lower_bounds, upper_bounds, strict=True)),
241 method="highs",
242 )
243 if not result.success:
244 msg = f"Could not find a maximum-return vertex (linear program: {result.message})"
245 raise ValueError(msg)
246 return np.asarray(result.x, dtype=np.float64)
249def _reject_degenerate_vertex(
250 a: NDArray[np.float64],
251 g: NDArray[np.float64],
252 free: NDArray[np.bool_],
253 active_ineq: NDArray[np.bool_],
254) -> None:
255 """Decline a maximum-return vertex whose free set cannot span the active rows.
257 The free set must span the equality rows together with the active inequality
258 rows: ``C = [A ; G_active]`` restricted to the free assets must have full row
259 rank, or the reduced KKT solve is singular. A degenerate maximum-return
260 vertex (a basic asset pinned on a bound) violates this; decline it with an
261 actionable diagnosis instead of letting it surface as an opaque "Singular
262 matrix" error downstream.
264 Args:
265 a: Equality-constraint matrix (``m x n``).
266 g: Inequality-constraint matrix (``p x n``); empty ``(0, n)`` when none.
267 free: Boolean mask of the assets strictly inside their box bounds.
268 active_ineq: Boolean mask of the tight (active) inequality rows.
270 Raises:
271 ValueError: If the free set does not span the active equality and
272 inequality rows.
273 """
274 c = np.vstack([a, g[active_ineq]])
275 mc = c.shape[0]
276 n_free = int(np.count_nonzero(free))
277 # rank(C[:, free]) <= min(mc, n_free), so fewer free assets than active rows
278 # is degenerate by itself. Testing this first also keeps matrix_rank off a
279 # zero-column block, whose empty singular-value reduction raises on numpy 2.0.
280 if n_free < mc or np.linalg.matrix_rank(c[:, free]) < mc:
281 msg = (
282 f"The maximum-return vertex is degenerate (free-set size {n_free}, "
283 f"active constraints {mc}): a basic asset sits exactly on a box bound, so the free set "
284 "does not span the active equality and inequality rows and the reduced KKT system is "
285 "singular. Tracing a frontier from a degenerate first vertex is not yet supported; perturb "
286 "the bounds or the constraints so the maximum-return vertex is non-degenerate."
287 )
288 raise ValueError(msg)
291def _free(
292 w: NDArray[np.float64], lower_bounds: NDArray[np.float64], upper_bounds: NDArray[np.float64]
293) -> NDArray[np.bool_]:
294 """Determine which asset should be free in the turning point.
296 This helper function identifies the asset that should be marked as free
297 in the turning point. It selects the asset that is furthest from its bounds,
298 which helps ensure numerical stability in the algorithm.
300 Args:
301 w: Vector of portfolio weights.
302 lower_bounds: Vector of lower bounds for asset weights.
303 upper_bounds: Vector of upper bounds for asset weights.
305 Returns:
306 A boolean vector indicating which asset is free (True) and which are blocked (False).
308 """
309 # Calculate the distance from each weight to its nearest bound
310 distance = np.min(np.array([np.abs(w - lower_bounds), np.abs(upper_bounds - w)]), axis=0)
312 # Find the index of the asset furthest from its bounds
313 index = np.argmax(distance)
315 # Create a boolean vector with only that asset marked as free
316 free = np.full_like(w, False, dtype=np.bool_)
317 free[index] = True
318 return free