Coverage for src/cvxcla/lasso.py: 100%
138 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-11 09:58 +0000
1"""LASSO / LARS regularisation path as a parametric active-set problem.
3This module shows that the Critical Line Algorithm's machinery is not specific to
4portfolios: the *same* ``cvxcla.pathtracer.trace`` loop, the *same*
5``QuadraticForm`` operator, and the *same* Bland event selection trace the LASSO
6homotopy. Only the problem-specific glue (the segment solve and what an event
7means) differs.
9The LASSO solves, for a response ``y`` and design matrix ``X``,
11 minimize 1/2 ||y - X beta||^2 + lam ||beta||_1
13and its minimiser ``beta(lam)`` is continuous and piecewise linear in the penalty
14``lam``. On a segment where the active set ``A`` (the support) and the signs
15``s_A`` are fixed,
17 beta_A(lam) = (X_A^T X_A)^{-1} (X_A^T y - lam s_A) = alpha_A - lam * beta_slope_A
18 correlation(lam) = X^T (y - X beta(lam)) = p + lam * q
20with ``|correlation_j| <= lam`` off the support and ``correlation_j = lam s_j`` on
21it. The role played by the covariance ``Sigma`` and mean ``mu`` in the CLA is
22played here by the Gram matrix ``H = X^T X`` (wrapped in ``DenseCovariance``) and
23the vector ``X^T y``.
25**Constraints.** Like the CLA, the path tracer admits general linear inequality
26constraints ``G beta <= h``. An active row enters the reduced KKT system exactly as
27in the CLA (the bordered Schur complement of ``cla.py``), and the generalised
28correlation that drives the enter/leave events carries the active-row multipliers,
29``correlation(lam) = X^T y - H beta(lam) - G_S^T eta(lam)``. The constrained path is
30still piecewise linear (a quadratic loss under a polyhedral penalty *and* polyhedral
31constraints; cf. Rosset and Zhu). We require ``h > 0`` so the path can start from
32``beta = 0`` with every row slack -- the same first vertex as the unconstrained
33LASSO. (Equality constraints, or ``h`` with a zero entry, need a feasibility seed
34analogous to the CLA's linear-programming first vertex, and are left to future work.)
36Event families, mirroring the CLA's "move to / leave a bound":
38* **leave** -- an active coefficient reaches zero: ``lam = alpha_j / beta_slope_j``.
39* **enter** -- an inactive (generalised) correlation reaches ``+/-lam``.
40* **activate** -- a slack inequality row's residual ``G_r beta - h_r`` reaches zero.
41* **release** -- an active row's multiplier ``eta_r`` reaches zero.
42"""
44from __future__ import annotations
46from dataclasses import dataclass, field
47from functools import cached_property
48from itertools import pairwise
49from typing import cast
51import numpy as np
52from numpy.typing import NDArray
54from ._lasso import LassoSegment, LassoState, scan_events, solve_segment
55from ._lasso_validate import validate_constraints, validate_design_inputs, validate_operator_inputs
56from .operators import DenseCovariance, GramCovariance, QuadraticForm
57from .pathtracer import InequalityConstrained, trace
60@dataclass(frozen=True)
61class Breakpoint:
62 """A vertex of the piecewise-linear LASSO path.
64 Attributes:
65 lam: The penalty value at this breakpoint.
66 beta: The coefficient vector ``beta(lam)``.
67 active: Boolean mask of the support (non-zero coefficients) on the
68 segment leaving this breakpoint towards smaller ``lam``.
69 """
71 lam: float
72 beta: NDArray[np.float64]
73 active: NDArray[np.bool_]
76@dataclass
77class Lasso(InequalityConstrained):
78 """The LASSO regularisation path, traced as a parametric active-set problem.
80 Constructing a ``Lasso`` traces the entire path from ``lam_max`` (where
81 ``beta = 0``) down to ``lam = 0`` (the least-squares fit on the final support,
82 subject to any active constraints), storing the breakpoints in ``path``. The
83 walk is driven by the same ``cvxcla.pathtracer.trace`` loop as the Critical Line
84 Algorithm.
86 Optional linear inequality constraints ``G beta <= h`` (with ``h > 0``) are
87 traced through the same bordered solve as the CLA's ``G w <= h`` rows.
89 The quadratic form may be given either as a dense design ``(x, y)`` (the usual
90 case, ``H = X^T X``) or, via :meth:`from_operator`, as a ``QuadraticForm``
91 operator with the linear term ``X^T y``. The operator route lets a structured
92 form, a diagonal-plus-low-rank factor model or a kernel, drive the path in
93 ``O(nk)`` per step without forming the ``n x n`` Gram matrix, exactly as on the
94 portfolio side.
96 Attributes:
97 x: Design matrix of shape ``(m, n)`` (``None`` in operator mode).
98 y: Response vector of shape ``(m,)`` (``None`` in operator mode).
99 g: Optional inequality matrix ``(p, n)`` of ``G beta <= h``; ``None`` means
100 the plain LASSO.
101 h: Optional inequality right-hand side ``(p,)``; must be strictly positive.
102 nonneg: When ``True``, restrict to the non-negative LASSO ``beta >= 0``;
103 the default ``False`` traces the ordinary signed path.
104 gram: When ``True``, drive the path with the ``GramCovariance`` data-matrix
105 backend (Woodbury solves in the ``m``-dimensional observation space),
106 never materialising the ``n x n`` Gram ``X^T X`` — the win in the
107 ``n >> m`` regime. The default ``False`` forms the dense Gram.
108 tol: Tolerance for event selection and the validity window.
109 path: The discovered breakpoints, populated on construction.
110 quad_form: Optional ``QuadraticForm`` operator ``H`` (operator mode).
111 linear: Optional linear term ``X^T y`` of shape ``(n,)`` (operator mode).
112 """
114 x: NDArray[np.float64] | None = None
115 y: NDArray[np.float64] | None = None
116 g: NDArray[np.float64] | None = None
117 h: NDArray[np.float64] | None = None
118 nonneg: bool = False # pragma: no mutate
119 gram: bool = False # pragma: no mutate
120 tol: float = 1e-9 # pragma: no mutate
121 path: list[Breakpoint] = field(default_factory=list)
122 quad_form: QuadraticForm | None = None # pragma: no mutate
123 linear: NDArray[np.float64] | None = None # pragma: no mutate
125 def __post_init__(self) -> None:
126 """Validate shapes and trace the full LASSO path.
128 Raises:
129 ValueError: If ``x`` is not 2d, ``y``'s length does not match ``x``, the
130 constraint shapes are inconsistent, or any ``h`` entry is not
131 strictly positive (which would make ``beta = 0`` infeasible).
132 """
133 if self.quad_form is not None or self.linear is not None:
134 self.linear = validate_operator_inputs(self.quad_form, self.linear, self.x, self.y)
135 else:
136 validate_design_inputs(self.x, self.y)
137 validate_constraints(self.g, self.h, self.dimension, self.tol)
138 trace(self)
140 @classmethod
141 def problem(cls, x: NDArray[np.float64], y: NDArray[np.float64]) -> LassoBuilder:
142 """Start a fluent :class:`cvxcla.builder.LassoBuilder` for a LASSO path.
144 The LASSO counterpart of :meth:`cvxcla.cla.CLA.problem`: chain
145 ``.inequality(G, h)`` and finish with ``.trace()``. The builder maps onto the
146 constructor arguments and adds no modelling power.
148 Args:
149 x: Design matrix of shape ``(m, n)``.
150 y: Response vector of shape ``(m,)``.
152 Returns:
153 A :class:`cvxcla.builder.LassoBuilder`.
154 """
155 return LassoBuilder(x, y)
157 @classmethod
158 def from_operator(
159 cls,
160 quad: QuadraticForm,
161 xty: NDArray[np.float64],
162 *,
163 g: NDArray[np.float64] | None = None,
164 h: NDArray[np.float64] | None = None,
165 nonneg: bool = False,
166 tol: float = 1e-9,
167 ) -> Lasso:
168 """Trace a LASSO path with the quadratic form supplied as an operator.
170 The regression counterpart of :class:`cvxcla.cla.CLA` accepting a
171 ``QuadraticForm`` covariance. Instead of a dense design ``X``, pass the Gram
172 operator ``H`` (anything implementing :class:`QuadraticForm`, for example a
173 :class:`cvxcla.operators.FactorCovariance` or a kernel) together with the
174 linear term ``X^T y``. The homotopy reaches ``H`` only through ``matvec`` and
175 ``solve_free``, so a diagonal-plus-low-rank factor model or a kernel traces
176 the path in ``O(nk)`` per step without ever forming the ``n x n`` matrix,
177 exactly as on the portfolio side. For the path to coincide with the
178 design-matrix LASSO one needs ``H = X^T X`` and ``xty = X^T y``
179 (Theorem 1); any positive-semidefinite operator whose free blocks are
180 positive definite traces a well-defined path.
182 Args:
183 quad: The quadratic form ``H`` as a :class:`QuadraticForm` operator.
184 xty: The linear term ``X^T y`` of shape ``(n,)``.
185 g: Optional inequality matrix of ``G beta <= h``.
186 h: Optional inequality right-hand side; entries must be strictly positive.
187 nonneg: Restrict the path to ``beta >= 0``.
188 tol: Tolerance for event selection and the validity window.
190 Returns:
191 A traced :class:`Lasso` whose ``path`` holds the breakpoints.
192 """
193 return cls(
194 quad_form=quad,
195 linear=np.asarray(xty, dtype=np.float64),
196 g=g,
197 h=h,
198 nonneg=nonneg,
199 tol=tol,
200 )
202 @cached_property
203 def quad(self) -> QuadraticForm:
204 """The Gram matrix ``X^T X`` as a ``QuadraticForm`` backend (cached: ``X`` is fixed).
206 With ``gram=True`` the data-matrix backend is used instead of forming the
207 ``n x n`` Gram: it solves through the Woodbury identity in the
208 ``m``-dimensional observation space and never materialises an ``n x n``
209 matrix, the win in the high-dimensional ``p >> n`` regime (more features than
210 observations). ``GramCovariance`` represents ``X_c^T X_c / (m-1)``, so scaling
211 the data by ``sqrt(m-1)`` recovers ``X^T X`` exactly for a **centred** design
212 (the standard LASSO convention; pass a column-centred ``x``).
213 """
214 if self.quad_form is not None:
215 return self.quad_form
216 # Not operator mode, so __post_init__ guarantees a design matrix.
217 x = cast("NDArray[np.float64]", self.x)
218 if self.gram:
219 m = x.shape[0]
220 return GramCovariance(x * np.sqrt(m - 1.0))
221 return DenseCovariance(x.T @ x)
223 @cached_property
224 def xty(self) -> NDArray[np.float64]:
225 """The linear data ``X^T y`` (the analogue of the CLA's expected returns; cached)."""
226 if self.linear is not None:
227 return self.linear
228 # Not operator mode, so __post_init__ guarantees a design (x, y).
229 x = cast("NDArray[np.float64]", self.x)
230 y = cast("NDArray[np.float64]", self.y)
231 return x.T @ y
233 @property
234 def dimension(self) -> int:
235 """Number of features ``n`` (the problem dimension for the path tracer)."""
236 if self.x is not None:
237 return int(self.x.shape[1])
238 return int(self.xty.shape[0])
240 @property
241 def lam_max(self) -> float:
242 """The smallest penalty at which ``beta = 0`` is optimal: ``||X^T y||_inf``.
244 With ``h > 0`` every inequality row is slack at ``beta = 0`` (zero
245 multiplier), so the unconstrained threshold is unchanged.
246 """
247 return float(np.max(np.abs(self.xty)))
249 def begin(self) -> tuple[float, LassoState]:
250 """Record the all-zero solution at the start penalty and enter the first coordinate.
252 For the plain or inequality-constrained LASSO the start is
253 ``lam_max = ||X^T y||_inf`` and the most-correlated coordinate enters with its
254 sign. Under the non-negative restriction ``beta >= 0`` the l1 penalty becomes
255 the linear term ``lam * 1^T beta``, only positive correlations can enter, so
256 the start is ``lam_max = max_j (X^T y)_j`` and the coordinate enters with sign
257 ``+``. When no coordinate can enter (e.g. every correlation is non-positive
258 under ``beta >= 0``), ``beta = 0`` is optimal for all ``lambda`` and the path
259 is the single point.
260 """
261 n = self.dimension
262 xty = self.xty
263 rows_active = np.zeros(self.g_matrix.shape[0], dtype=bool)
264 if self.nonneg:
265 lam_max = float(np.max(xty)) if n else 0.0
266 j0, s0 = int(np.argmax(xty)), 1.0
267 else:
268 lam_max = self.lam_max
269 j0 = int(np.argmax(np.abs(xty)))
270 s0 = float(np.sign(xty[j0]))
272 self.path.append(Breakpoint(max(lam_max, 0.0), np.zeros(n), np.zeros(n, dtype=bool)))
273 active = np.zeros(n, dtype=bool)
274 signs = np.zeros(n)
275 if lam_max > self.tol:
276 active[j0] = True
277 signs[j0] = s0
278 return max(lam_max, 0.0), LassoState(active, signs, rows_active, max(lam_max, 0.0))
280 def segment(self, state: LassoState) -> LassoSegment:
281 """Solve the affine segment for the current support, signs, and active rows.
283 Delegates to :func:`cvxcla._lasso.solve_segment`; see there for the
284 bordered Schur solve that admits active inequality rows.
285 """
286 return solve_segment(self.quad, self.xty, self.g_matrix, self.h_vector, state)
288 def event_matrix(self, state: LassoState, segment: LassoSegment) -> NDArray[np.float64]:
289 """Return the ``(n + p, 4)`` matrix of candidate critical lambdas.
291 Delegates to :func:`cvxcla._lasso.scan_events`; see there for the
292 coordinate (leave/enter) and inequality-row (activate/release) events.
293 """
294 return scan_events(self.dimension, self.g_matrix, self.h_vector, self.tol, self.nonneg, state, segment)
296 def step(self, state: LassoState, segment: LassoSegment, sec: int, direction: int, lam: float) -> LassoState:
297 """Record the breakpoint at ``lam`` after flipping coordinate or row ``sec``.
299 For a coordinate (``sec < n``): direction 0 removes it from the support, 1/2
300 add it with sign ``+1``/``-1``. For an inequality row (``sec >= n``):
301 direction 0 activates the row, 1 releases it. The path is continuous across
302 the flip, so the recorded coefficients are the old segment at ``lam``.
303 """
304 n = self.dimension
305 active = state.active.copy()
306 signs = state.signs.copy()
307 rows_active = state.rows_active.copy()
308 if sec < n:
309 if direction == 0:
310 active[sec] = False
311 signs[sec] = 0.0
312 else:
313 active[sec] = True
314 signs[sec] = 1.0 if direction == 1 else -1.0
315 else:
316 rows_active[sec - n] = direction == 0
318 beta = segment.alpha - lam * segment.beta_slope
319 self.path.append(Breakpoint(lam, beta, active.copy()))
320 return LassoState(active, signs, rows_active, lam)
322 def finish(self, state: LassoState, segment: LassoSegment) -> None:
323 """Record the ``lam = 0`` endpoint: the least-squares fit on the final support."""
324 self.path.append(Breakpoint(0.0, segment.alpha.copy(), state.active.copy()))
326 def solution(self, lam: float) -> NDArray[np.float64]:
327 """Evaluate the piecewise-linear path at penalty ``lam``.
329 Args:
330 lam: The penalty value at which to evaluate ``beta``.
332 Returns:
333 The coefficient vector ``beta(lam)``, by linear interpolation between
334 the bracketing breakpoints (clamped to the path's endpoints).
335 """
336 ordered = sorted(self.path, key=lambda bp: bp.lam)
337 if lam <= ordered[0].lam:
338 return ordered[0].beta
339 if lam >= ordered[-1].lam:
340 return ordered[-1].beta
341 for lo, hi in pairwise(ordered):
342 if lo.lam <= lam <= hi.lam:
343 weight = (lam - lo.lam) / (hi.lam - lo.lam)
344 return (1.0 - weight) * lo.beta + weight * hi.beta
345 msg = "lam lies within the path range but no bracketing segment was found" # pragma: no cover
346 raise AssertionError(msg) # pragma: no cover
349class LassoBuilder:
350 """Chainable builder for a LASSO regularisation-path problem.
352 The LASSO counterpart of :class:`cvxcla.cla.ProblemBuilder`. Construct one via
353 :meth:`Lasso.problem`, optionally add inequality constraints with
354 :meth:`inequality`, and finish with :meth:`trace`, which builds the
355 :class:`Lasso` and traces the entire regularisation path. Like the CLA builder
356 it adds no modelling power: it accepts the same ``G beta <= h`` rows the
357 ``Lasso`` already supports and nothing else.
359 Examples:
360 >>> import numpy as np
361 >>> from cvxcla import Lasso
362 >>> rng = np.random.default_rng(0)
363 >>> x = rng.standard_normal((30, 5))
364 >>> y = rng.standard_normal(30)
365 >>> lasso = Lasso.problem(x, y).trace()
366 >>> len(lasso.path) > 0
367 True
368 """
370 def __init__(self, x: NDArray[np.float64], y: NDArray[np.float64]) -> None:
371 """Start a builder for design matrix ``x`` and response ``y``.
373 Args:
374 x: Design matrix of shape ``(m, n)``.
375 y: Response vector of shape ``(m,)``.
376 """
377 self.x = np.asarray(x, dtype=np.float64)
378 self.y = np.asarray(y, dtype=np.float64)
379 self._g_blocks: list[NDArray[np.float64]] = []
380 self._h_blocks: list[NDArray[np.float64]] = []
381 self._nonneg = False
383 def non_negative(self) -> LassoBuilder:
384 """Restrict the coefficients to ``beta >= 0`` (the non-negative LASSO).
386 Under ``beta >= 0`` the l1 penalty collapses to the linear term
387 ``lam * sum(beta)``, so the path is the standard one restricted to positive
388 signs -- structurally the CLA's box-bounded parametric QP.
390 Returns:
391 ``self``, for chaining.
392 """
393 self._nonneg = True
394 return self
396 def inequality(self, g: NDArray[np.float64], h: float | NDArray[np.float64]) -> LassoBuilder:
397 """Add one or more inequality rows ``G beta <= h`` (repeated calls accumulate).
399 Args:
400 g: A length-``n`` row vector or a ``(p, n)`` matrix.
401 h: The matching right-hand side: a scalar for a single row, or a
402 length-``p`` vector. Each entry must be strictly positive (so
403 ``beta = 0`` stays feasible), checked when the path is traced.
405 Returns:
406 ``self``, for chaining.
408 Raises:
409 ValueError: If ``g``'s column count is not ``n`` or ``h``'s length does
410 not match the rows of ``g``.
411 """
412 g_block = np.atleast_2d(np.asarray(g, dtype=np.float64))
413 h_block = np.atleast_1d(np.asarray(h, dtype=np.float64))
414 if self.x.ndim == 2 and g_block.shape[1] != self.x.shape[1]:
415 msg = f"inequality: coefficient matrix must have {self.x.shape[1]} columns, got shape {g_block.shape}"
416 raise ValueError(msg)
417 if h_block.shape[0] != g_block.shape[0]:
418 msg = f"inequality: h must have {g_block.shape[0]} entries to match the rows, got {h_block.shape[0]}"
419 raise ValueError(msg)
420 self._g_blocks.append(g_block)
421 self._h_blocks.append(h_block)
422 return self
424 def trace(self) -> Lasso:
425 """Assemble the pieces, build the ``Lasso``, and trace the full path.
427 Returns:
428 The traced :class:`Lasso`, whose ``path`` holds the breakpoints of the
429 (constrained) regularisation path.
430 """
431 g = np.vstack(self._g_blocks) if self._g_blocks else None
432 h = np.concatenate(self._h_blocks) if self._h_blocks else None
433 return Lasso(x=self.x, y=self.y, g=g, h=h, nonneg=self._nonneg)