Coverage for src/cvxcla/_lasso.py: 100%
66 statements
« prev ^ index » next coverage.py v7.15.4, created at 2026-08-11 09:58 +0000
1"""Per-segment numeric kernels for the LASSO regularisation path.
3The LASSO homotopy is driven by the same generic ``cvxcla.pathtracer.trace`` loop
4as the Critical Line Algorithm; the two problem-specific numeric kernels live here
5as pure functions, mirroring the CLA's ``cvxcla._kkt`` (the segment solve) and
6``cvxcla._events`` (the critical-lambda scan). :func:`solve_segment` solves the
7affine path valid on one segment; :func:`scan_events` stacks every candidate
8critical lambda into the ``(n + p, 4)`` matrix the tracer scans.
10Both are pure: they take arrays (and the ``QuadraticForm`` operator) in and return
11arrays, so ``cvxcla.lasso.Lasso`` reduces to the thin ``ParametricProblem`` glue
12that wires them into the tracer. The two ``NamedTuple`` carriers passed between the
13tracer and these kernels -- :class:`LassoState` (the vertex) and
14:class:`LassoSegment` (the affine piece) -- are defined here alongside the kernels
15that produce and consume them.
16"""
18from __future__ import annotations
20from typing import NamedTuple
22import numpy as np
23from numpy.typing import NDArray
25from .operators import QuadraticForm, bordered_solve
28class LassoState(NamedTuple):
29 """The support, sign pattern, active inequality rows, and current penalty.
31 ``lam`` is the penalty at the segment's upper end. The event scan uses it to
32 require *strict* progress to a smaller penalty: right after a coefficient
33 enters it sits at zero, so its leave event lies at the current ``lam``; without
34 the strict window the shared selector would re-fire it and the walk would cycle
35 between entering and leaving the same coordinate.
36 """
38 active: NDArray[np.bool_]
39 signs: NDArray[np.float64]
40 rows_active: NDArray[np.bool_]
41 lam: float
44class LassoSegment(NamedTuple):
45 """The affine path ``beta(lam) = alpha - lam * beta_slope`` and its multipliers.
47 ``eta_alpha``/``eta_slope`` give the active-row multiplier path
48 ``eta(lam) = eta_alpha + lam * eta_slope`` (length ``p``, nonzero only on active
49 rows); ``p``/``q`` give the generalised correlation ``p + lam * q``.
50 """
52 alpha: NDArray[np.float64]
53 beta_slope: NDArray[np.float64]
54 p: NDArray[np.float64]
55 q: NDArray[np.float64]
56 eta_alpha: NDArray[np.float64]
57 eta_slope: NDArray[np.float64]
60def solve_segment(
61 quad: QuadraticForm,
62 xty: NDArray[np.float64],
63 g_matrix: NDArray[np.float64],
64 h_vector: NDArray[np.float64],
65 state: LassoState,
66) -> LassoSegment:
67 """Solve the affine segment for the current support, signs, and active rows.
69 With no active rows this is the plain LASSO solve against the Gram
70 submatrix. With active rows it is the bordered Schur solve of the CLA: the
71 active rows ``G_S`` enter the reduced KKT system as extra equality rows.
73 Args:
74 quad: The quadratic form ``H`` (``X^T X``) as a :class:`QuadraticForm`.
75 xty: The linear term ``X^T y`` of shape ``(n,)``.
76 g_matrix: Inequality matrix ``G`` of ``G beta <= h`` (``(p, n)``).
77 h_vector: Inequality right-hand side ``h`` (length ``p``).
78 state: The current support, signs, and active-row masks.
80 Returns:
81 The :class:`LassoSegment` affine path and its multipliers.
82 """
83 n = xty.shape[0]
84 active, signs, rows_active = state.active, state.signs, state.rows_active
85 alpha = np.zeros(n)
86 beta_slope = np.zeros(n)
87 eta_alpha = np.zeros(g_matrix.shape[0])
88 eta_slope = np.zeros(g_matrix.shape[0])
90 xty_s = xty[active]
91 signs_s = signs[active]
92 if not np.any(active):
93 # Empty support (e.g. the non-negative path when no correlation is
94 # positive): beta = 0, correlation = X^T y, and there is nothing to solve.
95 return LassoSegment(alpha, beta_slope, xty.copy(), np.zeros(n), eta_alpha, eta_slope)
97 # The active rows G_RS act as equality rows in the reduced KKT system, exactly
98 # the CLA's bordered Schur solve (operators.bordered_solve). With no active rows
99 # (|R| = 0) this reduces to the plain LASSO solve beta_S(lam) = H_SS^{-1}(xty_S -
100 # lam s_S). The slope's constraint right-hand side is zero, and the slope
101 # multiplier nu carries the opposite sign convention to eta(lam) (beta = alpha -
102 # lam beta_slope here, vs w = r_alpha + lam r_beta in the CLA), hence the flip.
103 g_rs = g_matrix[np.ix_(rows_active, active)] # |R| x |S|
104 h_r = h_vector[rows_active]
105 x_const, x_slope, eta_a, nu_slope = bordered_solve(quad, active, g_rs, xty_s, signs_s, h_r, np.zeros(g_rs.shape[0]))
106 alpha[active] = x_const
107 beta_slope[active] = x_slope
108 eta_alpha[rows_active] = eta_a
109 eta_slope[rows_active] = -nu_slope
111 # Generalised correlation c(lam) = xty - H beta(lam) - G_R^T eta(lam) = p + lam q.
112 g_r = g_matrix[rows_active]
113 eta_a_r = eta_alpha[rows_active]
114 eta_s_r = eta_slope[rows_active]
115 p = xty - quad.matvec(alpha) - g_r.T @ eta_a_r
116 q = quad.matvec(beta_slope) - g_r.T @ eta_s_r
117 return LassoSegment(alpha, beta_slope, p, q, eta_alpha, eta_slope)
120def scan_events(
121 n: int,
122 g_matrix: NDArray[np.float64],
123 h_vector: NDArray[np.float64],
124 tol: float,
125 nonneg: bool,
126 state: LassoState,
127 segment: LassoSegment,
128) -> NDArray[np.float64]:
129 """Return the ``(n + p, 4)`` matrix of candidate critical lambdas.
131 Rows ``0..n-1`` are coordinate events (col 0 leave, col 1 enter ``+``, col 2
132 enter ``-``); rows ``n..n+p-1`` are inequality-row events (col 0 activate, col
133 1 release). Entries are ``-inf`` where the event cannot occur, and only
134 events strictly inside ``(tol, lam - tol)`` are kept.
136 Args:
137 n: The problem dimension (number of features).
138 g_matrix: Inequality matrix ``G`` of ``G beta <= h`` (``(p, n)``).
139 h_vector: Inequality right-hand side ``h`` (length ``p``).
140 tol: Tolerance for event selection and the validity window.
141 nonneg: When ``True`` the enter-negative events are disabled (``beta >= 0``).
142 state: The current support, signs, and active-row masks.
143 segment: The affine path and multipliers from :func:`solve_segment`.
145 Returns:
146 The ``(n + p, 4)`` matrix of critical lambdas.
147 """
148 rows = g_matrix.shape[0]
149 active = state.active
150 inactive = ~active
151 alpha, beta_slope, p, q, eta_alpha, eta_slope = segment
152 rows_active = state.rows_active
154 l_mat = np.full((n + rows, 4), -np.inf)
156 # leave: alpha_j - lam beta_slope_j = 0
157 leaves = active & (np.abs(beta_slope) > tol) # pragma: no mutate
158 l_mat[:n][leaves, 0] = alpha[leaves] / beta_slope[leaves]
160 # enter (+): p_j + lam q_j = +lam -> lam = p_j / (1 - q_j)
161 denom_pos = 1.0 - q
162 enters_pos = inactive & (np.abs(denom_pos) > tol) # pragma: no mutate
163 l_mat[:n][enters_pos, 1] = p[enters_pos] / denom_pos[enters_pos]
165 # enter (-): p_j + lam q_j = -lam -> lam = -p_j / (1 + q_j). Disabled under the
166 # non-negative restriction beta >= 0, where only positive entries are allowed.
167 if not nonneg:
168 denom_neg = 1.0 + q
169 enters_neg = inactive & (np.abs(denom_neg) > tol) # pragma: no mutate
170 l_mat[:n][enters_neg, 2] = -p[enters_neg] / denom_neg[enters_neg]
172 if rows:
173 slope_row = g_matrix @ beta_slope # d/d(-lam) of the row value
174 level_row = g_matrix @ alpha - h_vector # G_r alpha - h_r
175 # activate: the row value G_r beta(lam) = level + h_r - lam slope rises to
176 # the cap h_r as lam decreases when its slope d(value)/d(-lam) = slope_row
177 # is positive; the crossing is lam = (G_r alpha - h_r) / (G_r beta_slope).
178 inactive_rows = ~rows_active & (slope_row > tol) # pragma: no mutate
179 l_mat[n:][inactive_rows, 0] = level_row[inactive_rows] / slope_row[inactive_rows]
180 # release: eta_r(lam) = eta_alpha + lam eta_slope -> 0 from eta > 0, i.e. eta_slope > 0.
181 releasing = rows_active & (eta_slope > tol) # pragma: no mutate
182 l_mat[n:][releasing, 1] = -eta_alpha[releasing] / eta_slope[releasing]
184 # Keep only events that make strict progress to a smaller, positive penalty.
185 l_mat[(l_mat <= tol) | (l_mat >= state.lam - tol)] = -np.inf
186 return l_mat