Coverage for src/cvxcla/lasso.py: 100%

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 two are closer than analogues. Under ``Sigma = X^T X`` and ``mu = X^T y`` the

10constrained LASSO path ``beta(lam)`` solves the gross-exposure-capped Markowitz

11program ``min 1/2 w^T Sigma w - mu^T w`` s.t. ``||w||_1 <= c`` at

12``c = ||beta(lam)||_1``, under the same linear constraints, and the two paths share

13their breakpoints wherever ``c`` is strictly decreasing (Schmelzer and Hastie,

14"The Critical Line Algorithm and the Constrained LASSO: One Curve, Two

15Literatures", arXiv:2609.25704, Theorem 1). With homogeneous constraints, the

16tilt sweep of ``CLA(leverage=c)`` is that path rescaled,

17``w_c(lam) = lam * beta`` with ``||beta||_1 = c / lam`` (Corollary 2).

19The LASSO solves, for a response ``y`` and design matrix ``X``,

21 minimize 1/2 ||y - X beta||^2 + lam ||beta||_1

23and its minimiser ``beta(lam)`` is continuous and piecewise linear in the penalty

24``lam``. On a segment where the active set ``A`` (the support) and the signs

25``s_A`` are fixed,

27 beta_A(lam) = (X_A^T X_A)^{-1} (X_A^T y - lam s_A) = alpha_A - lam * beta_slope_A

28 correlation(lam) = X^T (y - X beta(lam)) = p + lam * q

30with ``|correlation_j| <= lam`` off the support and ``correlation_j = lam s_j`` on

31it. The role played by the covariance ``Sigma`` and mean ``mu`` in the CLA is

32played here by the Gram matrix ``H = X^T X`` (wrapped in ``DenseCovariance``) and

33the vector ``X^T y``.

35**Constraints.** Like the CLA, the path tracer admits general linear inequality

36constraints ``G beta <= h``. An active row enters the reduced KKT system exactly as

37in the CLA (the bordered Schur complement of ``cla.py``), and the generalised

38correlation that drives the enter/leave events carries the active-row multipliers,

39``correlation(lam) = X^T y - H beta(lam) - G_S^T eta(lam)``. The constrained path is

40still piecewise linear (a quadratic loss under a polyhedral penalty *and* polyhedral

41constraints; cf. Rosset and Zhu). We require ``h > 0`` so the path can start from

42``beta = 0`` with every row slack -- the same first vertex as the unconstrained

43LASSO. Homogeneous equality constraints ``A beta = 0`` are traced by a different

44route: one leverage-capped CLA, rescaled (see :mod:`cvxcla._lasso_cla`).

46Event families, mirroring the CLA's "move to / leave a bound":

48* **leave** -- an active coefficient reaches zero: ``lam = alpha_j / beta_slope_j``.

49* **enter** -- an inactive (generalised) correlation reaches ``+/-lam``.

50* **activate** -- a slack inequality row's residual ``G_r beta - h_r`` reaches zero.

51* **release** -- an active row's multiplier ``eta_r`` reaches zero.

52"""

54from __future__ import annotations

56from bisect import bisect_right

57from dataclasses import dataclass, field

58from functools import cached_property

59from typing import cast

61import numpy as np

62from numpy.typing import NDArray

64from ._builders import LassoBuilder

65from ._lasso import LassoSegment, LassoState, scan_events, solve_segment

66from ._lasso_cla import equality_path

67from ._lasso_validate import (

68 validate_constraints,

69 validate_design_inputs,

70 validate_equality,

71 validate_operator_inputs,

72)

73from .operators import DenseCovariance, GramCovariance, QuadraticForm

74from .pathtracer import InequalityConstrained, trace

77@dataclass(frozen=True)

78class Breakpoint:

79 """A vertex of the piecewise-linear LASSO path.

81 Attributes:

82 lam: The penalty value at this breakpoint.

83 beta: The coefficient vector ``beta(lam)``.

84 active: Boolean mask of the support (non-zero coefficients) on the

85 segment leaving this breakpoint towards smaller ``lam``.

86 """

88 lam: float

89 beta: NDArray[np.float64]

90 active: NDArray[np.bool_]

93@dataclass

94class Lasso(InequalityConstrained):

95 """The LASSO regularisation path, traced as a parametric active-set problem.

97 Constructing a ``Lasso`` traces the entire path from ``lam_max`` (where

98 ``beta = 0``) down to ``lam = 0`` (the least-squares fit on the final support,

99 subject to any active constraints), storing the breakpoints in ``path``. The

100 walk is driven by the same ``cvxcla.pathtracer.trace`` loop as the Critical Line

101 Algorithm.

102

103 Optional linear inequality constraints ``G beta <= h`` (with ``h > 0``) are

104 traced through the same bordered solve as the CLA's ``G w <= h`` rows.

105 Homogeneous equality constraints ``A beta = 0`` (sum-to-zero, contrasts) are

106 traced instead through one leverage-capped CLA under ``Sigma = X^T X`` and

107 ``mu = X^T y``, which traces the same curve (Schmelzer and Hastie,

108 arXiv:2609.25704, Theorem 1 and Corollary 2); see :mod:`cvxcla._lasso_cla`.

109

110 The quadratic form may be given either as a dense design ``(x, y)`` (the usual

111 case, ``H = X^T X``) or, via :meth:`from_operator`, as a ``QuadraticForm``

112 operator with the linear term ``X^T y``. The operator route lets a structured

113 form, a diagonal-plus-low-rank factor model or a kernel, drive the path in

114 ``O(nk)`` per step without forming the ``n x n`` Gram matrix, exactly as on the

115 portfolio side.

116

117 Attributes:

118 x: Design matrix of shape ``(m, n)`` (``None`` in operator mode).

119 y: Response vector of shape ``(m,)`` (``None`` in operator mode).

120 g: Optional inequality matrix ``(p, n)`` of ``G beta <= h``; ``None`` means

121 the plain LASSO.

122 h: Optional inequality right-hand side ``(p,)``; must be strictly positive.

123 nonneg: When ``True``, restrict to the non-negative LASSO ``beta >= 0``;

124 the default ``False`` traces the ordinary signed path.

125 gram: When ``True``, drive the path with the ``GramCovariance`` data-matrix

126 backend (Woodbury solves in the ``m``-dimensional observation space),

127 never materialising the ``n x n`` Gram ``X^T X`` — the win in the

128 ``n >> m`` regime. The default ``False`` forms the dense Gram.

129 tol: Tolerance for event selection and the validity window.

130 path: The discovered breakpoints, populated on construction.

131 quad_form: Optional ``QuadraticForm`` operator ``H`` (operator mode).

132 linear: Optional linear term ``X^T y`` of shape ``(n,)`` (operator mode).

133 a: Optional equality matrix ``(m, n)`` of ``A beta = 0``. It cannot be

134 combined with ``g``/``h``, and it needs a positive-definite Gram, so that

135 the constrained least-squares end of the path is unique.

136 """

137

138 x: NDArray[np.float64] | None = None

139 y: NDArray[np.float64] | None = None

140 g: NDArray[np.float64] | None = None

141 h: NDArray[np.float64] | None = None

142 nonneg: bool = False # pragma: no mutate

143 gram: bool = False # pragma: no mutate

144 tol: float = 1e-9 # pragma: no mutate

145 path: list[Breakpoint] = field(default_factory=list)

146 quad_form: QuadraticForm | None = None # pragma: no mutate

147 linear: NDArray[np.float64] | None = None # pragma: no mutate

148 a: NDArray[np.float64] | None = None # pragma: no mutate

149

150 def __post_init__(self) -> None:

151 """Validate shapes and trace the full LASSO path.

152

153 Raises:

154 ValueError: If ``x`` is not 2d, ``y``'s length does not match ``x``, the

155 constraint shapes are inconsistent, any ``h`` entry is not strictly

156 positive (which would make ``beta = 0`` infeasible), or equality

157 rows ``a`` are combined with ``g`` or cannot be traced.

158 """

159 if self.quad_form is not None or self.linear is not None:

160 self.linear = validate_operator_inputs(self.quad_form, self.linear, self.x, self.y)

161 else:

162 validate_design_inputs(self.x, self.y)

163 validate_constraints(self.g, self.h, self.dimension, self.tol)

164 if self.a is not None:

165 self.a = validate_equality(self.a, self.g, self.dimension)

166 path = equality_path(self.quad, self.xty, self.a, self.nonneg, self.tol)

167 self.path.extend(Breakpoint(lam, beta, active) for lam, beta, active in path)

168 return

169 trace(self)

170

171 @classmethod

172 def problem(cls, x: NDArray[np.float64], y: NDArray[np.float64]) -> LassoBuilder[Lasso]:

173 """Start a fluent :class:`cvxcla._builders.LassoBuilder` for a LASSO path.

174

175 The LASSO counterpart of :meth:`cvxcla.cla.CLA.problem`: chain

176 ``.inequality(G, h)`` or ``.equality(A)`` and finish with ``.trace()``. The builder maps onto the

177 constructor arguments and adds no modelling power.

178

179 ``cls`` is handed to the builder as the class it should construct, which is

180 what lets the builder live in a leaf module that never imports this one.

181

182 Args:

183 x: Design matrix of shape ``(m, n)``.

184 y: Response vector of shape ``(m,)``.

185

186 Returns:

187 A :class:`cvxcla._builders.LassoBuilder` whose ``.trace()`` returns a

188 traced ``Lasso``.

189 """

190 return LassoBuilder(x, y, solver=cls)

191

192 @classmethod

193 def from_operator(

194 cls,

195 quad: QuadraticForm,

196 xty: NDArray[np.float64],

197 *,

198 g: NDArray[np.float64] | None = None,

199 h: NDArray[np.float64] | None = None,

200 a: NDArray[np.float64] | None = None,

201 nonneg: bool = False,

202 tol: float = 1e-9,

203 ) -> Lasso:

204 """Trace a LASSO path with the quadratic form supplied as an operator.

205

206 The regression counterpart of :class:`cvxcla.cla.CLA` accepting a

207 ``QuadraticForm`` covariance. Instead of a dense design ``X``, pass the Gram

208 operator ``H`` (anything implementing :class:`QuadraticForm`, for example a

209 :class:`cvxcla.operators.FactorCovariance` or a kernel) together with the

210 linear term ``X^T y``. The homotopy reaches ``H`` only through ``matvec`` and

211 ``solve_free``, so a diagonal-plus-low-rank factor model or a kernel traces

212 the path in ``O(nk)`` per step without ever forming the ``n x n`` matrix,

213 exactly as on the portfolio side. For the path to coincide with the

214 design-matrix LASSO one needs ``H = X^T X`` and ``xty = X^T y``

215 (Theorem 1); any positive-semidefinite operator whose free blocks are

216 positive definite traces a well-defined path.

217

218 Args:

219 quad: The quadratic form ``H`` as a :class:`QuadraticForm` operator.

220 xty: The linear term ``X^T y`` of shape ``(n,)``.

221 g: Optional inequality matrix of ``G beta <= h``.

222 h: Optional inequality right-hand side; entries must be strictly positive.

223 a: Optional equality matrix of ``A beta = 0`` (not combined with ``g``).

224 nonneg: Restrict the path to ``beta >= 0``.

225 tol: Tolerance for event selection and the validity window.

226

227 Returns:

228 A traced :class:`Lasso` whose ``path`` holds the breakpoints.

229 """

230 return cls(

231 quad_form=quad,

232 linear=np.asarray(xty, dtype=np.float64),

233 g=g,

234 h=h,

235 a=a,

236 nonneg=nonneg,

237 tol=tol,

238 )

239

240 @cached_property

241 def quad(self) -> QuadraticForm:

242 """The Gram matrix ``X^T X`` as a ``QuadraticForm`` backend (cached: ``X`` is fixed).

243

244 With ``gram=True`` the data-matrix backend is used instead of forming the

245 ``n x n`` Gram: it solves through the Woodbury identity in the

246 ``m``-dimensional observation space and never materialises an ``n x n``

247 matrix, the win in the high-dimensional ``p >> n`` regime (more features than

248 observations). ``GramCovariance`` represents ``X_c^T X_c / (m-1)``, so scaling

249 the data by ``sqrt(m-1)`` recovers ``X^T X`` exactly for a **centred** design

250 (the standard LASSO convention; pass a column-centred ``x``).

251 """

252 if self.quad_form is not None:

253 return self.quad_form

254 # Not operator mode, so __post_init__ guarantees a design matrix.

255 x = cast("NDArray[np.float64]", self.x)

256 if self.gram:

257 m = x.shape[0]

258 return GramCovariance(x * np.sqrt(m - 1.0))

259 return DenseCovariance(x.T @ x)

260

261 @cached_property

262 def xty(self) -> NDArray[np.float64]:

263 """The linear data ``X^T y`` (the analogue of the CLA's expected returns; cached)."""

264 if self.linear is not None:

265 return self.linear

266 # Not operator mode, so __post_init__ guarantees a design (x, y).

267 x = cast("NDArray[np.float64]", self.x)

268 y = cast("NDArray[np.float64]", self.y)

269 return x.T @ y

270

271 @property

272 def dimension(self) -> int:

273 """Number of features ``n`` (the problem dimension for the path tracer)."""

274 if self.x is not None:

275 return int(self.x.shape[1])

276 return int(self.xty.shape[0])

277

278 @property

279 def lam_max(self) -> float:

280 """The smallest penalty at which ``beta = 0`` is optimal: ``||X^T y||_inf``.

281

282 With ``h > 0`` every inequality row is slack at ``beta = 0`` (zero

283 multiplier), so the unconstrained threshold is unchanged. Equality rows

284 ``A beta = 0`` do change it: the correlation that can enter is the part of

285 ``X^T y`` outside the row space of ``A``, so the threshold is read off the

286 traced path instead.

287 """

288 if self.a is not None:

289 return self.path[0].lam

290 return float(np.max(np.abs(self.xty)))

291

292 def begin(self) -> tuple[float, LassoState]:

293 """Record the all-zero solution at the start penalty and enter the first coordinate.

294

295 For the plain or inequality-constrained LASSO the start is

296 ``lam_max = ||X^T y||_inf`` and the most-correlated coordinate enters with its

297 sign. Under the non-negative restriction ``beta >= 0`` the l1 penalty becomes

298 the linear term ``lam * 1^T beta``, only positive correlations can enter, so

299 the start is ``lam_max = max_j (X^T y)_j`` and the coordinate enters with sign

300 ``+``. When no coordinate can enter (e.g. every correlation is non-positive

301 under ``beta >= 0``), ``beta = 0`` is optimal for all ``lambda`` and the path

302 is the single point.

303 """

304 n = self.dimension

305 xty = self.xty

306 rows_active = np.zeros(self.g_matrix.shape[0], dtype=bool)

307 if self.nonneg:

308 lam_max = float(np.max(xty)) if n else 0.0

309 j0, s0 = int(np.argmax(xty)), 1.0

310 else:

311 lam_max = self.lam_max

312 j0 = int(np.argmax(np.abs(xty)))

313 s0 = float(np.sign(xty[j0]))

314

315 self.path.append(Breakpoint(max(lam_max, 0.0), np.zeros(n), np.zeros(n, dtype=bool)))

316 active = np.zeros(n, dtype=bool)

317 signs = np.zeros(n)

318 if lam_max > self.tol:

319 active[j0] = True

320 signs[j0] = s0

321 return max(lam_max, 0.0), LassoState(active, signs, rows_active, max(lam_max, 0.0))

322

323 def segment(self, state: LassoState) -> LassoSegment:

324 """Solve the affine segment for the current support, signs, and active rows.

325

326 Delegates to :func:`cvxcla._lasso.solve_segment`; see there for the

327 bordered Schur solve that admits active inequality rows.

328 """

329 return solve_segment(self.quad, self.xty, self.g_matrix, self.h_vector, state)

330

331 def event_matrix(self, state: LassoState, segment: LassoSegment) -> NDArray[np.float64]:

332 """Return the ``(n + p, 4)`` matrix of candidate critical lambdas.

333

334 Delegates to :func:`cvxcla._lasso.scan_events`; see there for the

335 coordinate (leave/enter) and inequality-row (activate/release) events.

336 """

337 return scan_events(self.dimension, self.g_matrix, self.h_vector, self.tol, self.nonneg, state, segment)

338

339 def step(self, state: LassoState, segment: LassoSegment, sec: int, direction: int, lam: float) -> LassoState:

340 """Record the breakpoint at ``lam`` after flipping coordinate or row ``sec``.

341

342 For a coordinate (``sec < n``): direction 0 removes it from the support, 1/2

343 add it with sign ``+1``/``-1``. For an inequality row (``sec >= n``):

344 direction 0 activates the row, 1 releases it. The path is continuous across

345 the flip, so the recorded coefficients are the old segment at ``lam``.

346 """

347 n = self.dimension

348 active = state.active.copy()

349 signs = state.signs.copy()

350 rows_active = state.rows_active.copy()

351 if sec < n:

352 if direction == 0:

353 active[sec] = False

354 signs[sec] = 0.0

355 else:

356 active[sec] = True

357 signs[sec] = 1.0 if direction == 1 else -1.0

358 else:

359 rows_active[sec - n] = direction == 0

360

361 beta = segment.alpha - lam * segment.beta_slope

362 self.path.append(Breakpoint(lam, beta, active.copy()))

363 return LassoState(active, signs, rows_active, lam)

364

365 def finish(self, state: LassoState, segment: LassoSegment) -> None:

366 """Record the ``lam = 0`` endpoint: the least-squares fit on the final support."""

367 self.path.append(Breakpoint(0.0, segment.alpha.copy(), state.active.copy()))

368

369 def solution(self, lam: float) -> NDArray[np.float64]:

370 """Evaluate the piecewise-linear path at penalty ``lam``.

371

372 Args:

373 lam: The penalty value at which to evaluate ``beta``.

374

375 Returns:

376 The coefficient vector ``beta(lam)``, by linear interpolation between

377 the bracketing breakpoints (clamped to the path's endpoints).

378 """

379 ordered = sorted(self.path, key=lambda bp: bp.lam)

380 if lam <= ordered[0].lam:

381 return ordered[0].beta

382 if lam >= ordered[-1].lam:

383 return ordered[-1].beta

384 # Strictly inside the range, so lo.lam <= lam < hi.lam and hi.lam > lo.lam.

385 idx = bisect_right([bp.lam for bp in ordered], lam)

386 lo, hi = ordered[idx - 1], ordered[idx]

387 weight = (lam - lo.lam) / (hi.lam - lo.lam)

388 return (1.0 - weight) * lo.beta + weight * hi.beta