Coverage for src/cvxcla/first.py: 100%
48 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-11 09:58 +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 free = (weights > lower_bounds + tol) & (weights < upper_bounds - tol)
128 active_ineq = (g @ weights >= h - tol) if g.shape[0] else np.zeros(0, dtype=bool)
130 _reject_degenerate_vertex(a, g, free, active_ineq)
131 return TurningPoint(free=free, weights=weights, active_ineq=active_ineq)
134def _solve_max_return_lp(
135 mean: NDArray[np.float64],
136 lower_bounds: NDArray[np.float64],
137 upper_bounds: NDArray[np.float64],
138 a: NDArray[np.float64],
139 b: NDArray[np.float64],
140 g: NDArray[np.float64],
141 h: NDArray[np.float64],
142) -> NDArray[np.float64]:
143 """Solve the maximum-return linear program and return its vertex weights.
145 ``maximize mean @ w`` (as ``minimize -mean @ w``) subject to ``A w = b``,
146 ``G w <= h`` and the box bounds, via HiGHS. The inequality rows are passed
147 only when ``g`` is non-empty.
149 Args:
150 mean: Vector of expected returns.
151 lower_bounds: Lower box bounds.
152 upper_bounds: Upper box bounds.
153 a: Equality-constraint matrix (``m x n``).
154 b: Equality-constraint right-hand side (length ``m``).
155 g: Inequality-constraint matrix (``p x n``); empty ``(0, n)`` when none.
156 h: Inequality-constraint right-hand side (length ``p``).
158 Returns:
159 The vertex weights ``w`` as a 1d ``float64`` array.
161 Raises:
162 ValueError: If the linear program is infeasible or unbounded.
163 """
164 has_ineq = g.shape[0] > 0
165 result = linprog(
166 c=-np.asarray(mean, dtype=np.float64),
167 A_eq=np.asarray(a, dtype=np.float64),
168 b_eq=np.asarray(b, dtype=np.float64),
169 A_ub=g if has_ineq else None,
170 b_ub=h if has_ineq else None,
171 bounds=list(zip(lower_bounds, upper_bounds, strict=True)),
172 method="highs",
173 )
174 if not result.success:
175 msg = f"Could not find a maximum-return vertex (linear program: {result.message})"
176 raise ValueError(msg)
177 return np.asarray(result.x, dtype=np.float64)
180def _reject_degenerate_vertex(
181 a: NDArray[np.float64],
182 g: NDArray[np.float64],
183 free: NDArray[np.bool_],
184 active_ineq: NDArray[np.bool_],
185) -> None:
186 """Decline a maximum-return vertex whose free set cannot span the active rows.
188 The free set must span the equality rows together with the active inequality
189 rows: ``C = [A ; G_active]`` restricted to the free assets must have full row
190 rank, or the reduced KKT solve is singular. A degenerate maximum-return
191 vertex (a basic asset pinned on a bound) violates this; decline it with an
192 actionable diagnosis instead of letting it surface as an opaque "Singular
193 matrix" error downstream.
195 Args:
196 a: Equality-constraint matrix (``m x n``).
197 g: Inequality-constraint matrix (``p x n``); empty ``(0, n)`` when none.
198 free: Boolean mask of the assets strictly inside their box bounds.
199 active_ineq: Boolean mask of the tight (active) inequality rows.
201 Raises:
202 ValueError: If the free set does not span the active equality and
203 inequality rows.
204 """
205 c = np.vstack([a, g[active_ineq]])
206 mc = c.shape[0]
207 n_free = int(np.count_nonzero(free))
208 # rank(C[:, free]) <= min(mc, n_free), so fewer free assets than active rows
209 # is degenerate by itself. Testing this first also keeps matrix_rank off a
210 # zero-column block, whose empty singular-value reduction raises on numpy 2.0.
211 if n_free < mc or np.linalg.matrix_rank(c[:, free]) < mc:
212 msg = (
213 f"The maximum-return vertex is degenerate (free-set size {n_free}, "
214 f"active constraints {mc}): a basic asset sits exactly on a box bound, so the free set "
215 "does not span the active equality and inequality rows and the reduced KKT system is "
216 "singular. Tracing a frontier from a degenerate first vertex is not yet supported; perturb "
217 "the bounds or the constraints so the maximum-return vertex is non-degenerate."
218 )
219 raise ValueError(msg)
222def _free(
223 w: NDArray[np.float64], lower_bounds: NDArray[np.float64], upper_bounds: NDArray[np.float64]
224) -> NDArray[np.bool_]:
225 """Determine which asset should be free in the turning point.
227 This helper function identifies the asset that should be marked as free
228 in the turning point. It selects the asset that is furthest from its bounds,
229 which helps ensure numerical stability in the algorithm.
231 Args:
232 w: Vector of portfolio weights.
233 lower_bounds: Vector of lower bounds for asset weights.
234 upper_bounds: Vector of upper bounds for asset weights.
236 Returns:
237 A boolean vector indicating which asset is free (True) and which are blocked (False).
239 """
240 # Calculate the distance from each weight to its nearest bound
241 distance = np.min(np.array([np.abs(w - lower_bounds), np.abs(upper_bounds - w)]), axis=0)
243 # Find the index of the asset furthest from its bounds
244 index = np.argmax(distance)
246 # Create a boolean vector with only that asset marked as free
247 free = np.full_like(w, False, dtype=np.bool_)
248 free[index] = True
249 return free