Coverage for src/cvxcla/_kkt.py: 100%
53 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-09-29 05:24 +0000
1"""Reduced KKT machinery for a Critical Line Algorithm turning point.
3At each turning point the active set is identified (which box bounds are held,
4which assets are free) and the reduced KKT system is solved by block elimination
5to produce the affine critical-line segment ``w(lam) = r_alpha + lam * r_beta``.
6Both steps are pure functions of the problem data and the active set, so they
7live here rather than on the ``CLA`` class; the covariance only ever enters
8through the ``QuadraticForm`` interface, so structured backends never
9materialise an ``n x n`` matrix.
10"""
12from __future__ import annotations
14from typing import NamedTuple
16import numpy as np
17from numpy.typing import NDArray
19from .operators import QuadraticForm, bordered_solve, cross
20from .types import TurningPoint
23def active_set(
24 free: NDArray[np.bool_],
25 weights: NDArray[np.float64],
26 lower: NDArray[np.float64],
27 upper: NDArray[np.float64],
28 tol: float,
29) -> tuple[NDArray[np.bool_], NDArray[np.bool_], NDArray[np.bool_], NDArray[np.float64]]:
30 """Identify the active set at a turning point and the weights pinned to bounds.
32 A blocked asset sitting (to tolerance) on a bound is held fixed there and
33 excluded from the reduced KKT solve; every other asset is *in*. Returns the
34 upper-bound mask, the lower-bound mask, the in-set mask, and the full-length
35 vector of weights fixed at their bounds.
37 Args:
38 free: Boolean mask of the assets free at the turning point.
39 weights: The turning point's weight vector.
40 lower: Per-asset lower bounds.
41 upper: Per-asset upper bounds.
42 tol: Tolerance for classifying a weight as sitting on a bound.
44 Returns:
45 ``(at_upper, at_lower, free_in, fixed_weights)``.
47 Raises:
48 RuntimeError: If every asset is blocked, which makes the reduced system
49 singular.
50 """
51 blocked = ~free
52 if np.all(blocked):
53 msg = "All variables cannot be blocked"
54 raise RuntimeError(msg)
56 at_upper = blocked & (np.abs(weights - upper) <= tol) # pragma: no mutate
57 at_lower = blocked & (np.abs(weights - lower) <= tol) # pragma: no mutate
58 free_in = ~(at_upper | at_lower)
60 fixed_weights = np.zeros(len(weights))
61 fixed_weights[at_upper] = upper[at_upper]
62 fixed_weights[at_lower] = lower[at_lower]
63 return at_upper, at_lower, free_in, fixed_weights
66def solve_kkt(
67 cov: QuadraticForm,
68 mean: NDArray[np.float64],
69 a: NDArray[np.float64],
70 b: NDArray[np.float64],
71 g: NDArray[np.float64],
72 h: NDArray[np.float64],
73 free_in: NDArray[np.bool_],
74 fixed_weights: NDArray[np.float64],
75 active_ineq: NDArray[np.bool_],
76) -> tuple[
77 NDArray[np.float64],
78 NDArray[np.float64],
79 NDArray[np.float64],
80 NDArray[np.float64],
81 NDArray[np.float64],
82 NDArray[np.float64],
83]:
84 """Solve the reduced KKT system for the current critical-line segment.
86 Block elimination over the *stacked* constraint matrix ``C = [A ; G_S]``,
87 where ``G_S`` are the currently-active inequality rows held at equality.
88 Because an active inequality row enters the stationarity/feasibility system
89 exactly as an equality row does, the same elimination handles both: a
90 multi-right-hand-side solve against the free covariance block ``Sigma_FF``
91 (via the backend, so structured covariances never materialise an
92 ``n x n`` matrix) feeds an ``(m + |S|) x (m + |S|)`` Schur complement
93 ``C_F Sigma_FF^{-1} C_F.T``. With no active inequality rows this is the
94 plain equality solve.
96 Args:
97 cov: The covariance as a ``QuadraticForm`` backend.
98 mean: Vector of expected returns.
99 a: Equality-constraint matrix ``A`` of ``A w = b``.
100 b: Equality-constraint right-hand side ``b``.
101 g: Inequality-constraint matrix ``G`` of ``G w <= h`` (``(p, n)``).
102 h: Inequality-constraint right-hand side ``h`` (length ``p``).
103 free_in: Boolean mask of the assets in the reduced solve.
104 fixed_weights: Full-length weights of the assets held at their bounds.
105 active_ineq: Boolean mask (length ``p``) of the active inequality rows.
107 Returns:
108 ``(r_alpha, r_beta, gamma, delta, eta_alpha, eta_beta)``: the affine
109 segment ``w(lam) = r_alpha + lam * r_beta``, the box-multiplier
110 gradients ``gamma``/``delta`` that drive the leave-a-bound events, and
111 the affine inequality multipliers ``eta_alpha + lam * eta_beta``
112 (length ``p``, non-zero only on active rows) that drive the
113 release-a-row events.
114 """
115 m = a.shape[0]
116 ns = len(mean)
117 p = g.shape[0]
118 out = ~free_in
119 # Stack the active inequality rows beneath the equality rows; the active
120 # rows are held at equality (g_i w = h_i), so C/d is the equality system
121 # of the reduced QP at this vertex.
122 c = np.vstack([a, g[active_ineq]])
123 d = np.concatenate([b, h[active_ineq]])
124 c_free = c[:, free_in]
126 # The reduced KKT system is the shared bordered solve: the constant system
127 # carries the blocked-weight shift -Sigma_FB w_B and the reduced constraint
128 # right-hand side d - C_B w_B; the slope system carries the mean with a zero
129 # constraint right-hand side (A r_beta = 0).
130 x_alpha, x_beta, nu_alpha, nu_beta = bordered_solve(
131 cov,
132 free_in,
133 c_free,
134 -cross(cov, free_in, fixed_weights),
135 mean[free_in],
136 d - c[:, out] @ fixed_weights[out],
137 np.zeros(c.shape[0]),
138 )
140 r_alpha = fixed_weights.copy()
141 r_alpha[free_in] = x_alpha
142 r_beta = np.zeros(ns)
143 r_beta[free_in] = x_beta
145 gamma = cov.matvec(r_alpha) + c.T @ nu_alpha
146 delta = cov.matvec(r_beta) + c.T @ nu_beta - mean
148 # The tail of the stacked multiplier is the inequality multiplier eta(lam)
149 # = eta_alpha + lam eta_beta, scattered back to full length p (zero on the
150 # inactive rows, which have no release event).
151 eta_alpha = np.zeros(p)
152 eta_beta = np.zeros(p)
153 eta_alpha[active_ineq] = nu_alpha[m:]
154 eta_beta[active_ineq] = nu_beta[m:]
155 return r_alpha, r_beta, gamma, delta, eta_alpha, eta_beta
158class Segment(NamedTuple):
159 """The affine critical-line segment valid at one turning point.
161 Bundles the affine path ``w(lam) = r_alpha + lam * r_beta``, the multiplier
162 gradients ``gamma``/``delta`` that drive the leave-a-bound events, and the
163 active-set masks the event scan needs. This is what ``CLA.segment`` returns
164 to the generic path tracer.
166 For general inequality constraints ``G w <= h`` the segment also carries the
167 affine inequality multipliers ``eta_alpha + lam * eta_beta`` (one entry per
168 inequality row; meaningful for *active* rows, which release when the
169 multiplier crosses zero) and the active-row mask ``active_ineq``. The slacks
170 that drive an *inactive* row becoming active are recomputed from
171 ``r_alpha``/``r_beta`` directly in :func:`cvxcla._events.ineq_event_ratios`.
172 """
174 r_alpha: NDArray[np.float64]
175 r_beta: NDArray[np.float64]
176 gamma: NDArray[np.float64]
177 delta: NDArray[np.float64]
178 at_upper: NDArray[np.bool_]
179 at_lower: NDArray[np.bool_]
180 free_in: NDArray[np.bool_]
181 active_ineq: NDArray[np.bool_]
182 eta_alpha: NDArray[np.float64]
183 eta_beta: NDArray[np.float64]
186def critical_segment(
187 cov: QuadraticForm,
188 mean: NDArray[np.float64],
189 a: NDArray[np.float64],
190 b: NDArray[np.float64],
191 g: NDArray[np.float64],
192 h: NDArray[np.float64],
193 lower: NDArray[np.float64],
194 upper: NDArray[np.float64],
195 tol: float,
196 state: TurningPoint,
197) -> Segment:
198 """Solve the reduced KKT system for the critical-line segment at ``state``.
200 Composes :func:`active_set` (which bounds are held at the turning point) with
201 :func:`solve_kkt` (the block-eliminated reduced solve) and bundles the result
202 with the masks the event scan needs.
204 Args:
205 cov: The covariance as a ``QuadraticForm`` backend.
206 mean: Vector of expected returns.
207 a: Equality-constraint matrix ``A`` of ``A w = b``.
208 b: Equality-constraint right-hand side ``b``.
209 g: Inequality-constraint matrix ``G`` of ``G w <= h`` (``(p, n)``).
210 h: Inequality-constraint right-hand side ``h`` (length ``p``).
211 lower: Per-asset lower bounds.
212 upper: Per-asset upper bounds.
213 tol: Tolerance for classifying a weight as sitting on a bound.
214 state: The turning point the segment starts from.
216 Returns:
217 The :class:`Segment` valid below ``state``.
218 """
219 at_upper, at_lower, free_in, fixed_weights = active_set(state.free, state.weights, lower, upper, tol)
220 r_alpha, r_beta, gamma, delta, eta_alpha, eta_beta = solve_kkt(
221 cov, mean, a, b, g, h, free_in, fixed_weights, state.active_ineq
222 )
223 return Segment(r_alpha, r_beta, gamma, delta, at_upper, at_lower, free_in, state.active_ineq, eta_alpha, eta_beta)