Coverage for src/cvxcla/operators/_core.py: 100%
36 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-06-27 06:45 +0000
1"""Shared contract and solver primitives for the covariance backends.
3This module defines the :class:`QuadraticForm` protocol that every backend
4implements, the deterministic symmetric reciprocal-condition helper they share,
5and :func:`bordered_solve`, the bordered KKT solve used by both
6``cvxcla.cla`` and ``cvxcla.lasso``. The concrete backends live in sibling
7modules (``dense``, ``factor``, ``gram``) and are re-exported from the package
8root ``cvxcla.operators``.
9"""
11from __future__ import annotations
13from typing import Protocol, cast, runtime_checkable
15import numpy as np
16from numpy.typing import NDArray
18# When the cheap LAPACK 1-norm condition *estimate* clears the singularity floor
19# by this margin, the decision is conclusive; an estimate landing within the
20# margin falls back to the exact symmetric rcond. The margin generously covers
21# the gap between the 1-norm estimate and the true 2-norm reciprocal condition
22# number, so the boolean answer is identical to the exact comparison.
23_RCOND_ESTIMATE_MARGIN = 1.0e3
26def _rcond_symmetric(block: NDArray[np.float64]) -> float:
27 """Reciprocal 2-norm condition number of a symmetric (PSD) ``block``.
29 Returns a value in ``[0, 1]``: ``1`` is perfectly conditioned and ``0`` is
30 numerically singular. The condition number is read off the symmetric
31 eigenvalues (``lambda_min / lambda_max``), so the result is deterministic
32 and independent of the BLAS/LAPACK build, unlike the residual of a solve
33 against a singular block. An empty block is treated as well conditioned.
34 """
35 if block.shape[0] == 0:
36 return 1.0
37 eigenvalues = np.linalg.eigvalsh(block)
38 lam_max = float(eigenvalues[-1])
39 if lam_max <= 0.0:
40 return 0.0
41 return max(float(eigenvalues[0]), 0.0) / lam_max
44@runtime_checkable
45class QuadraticForm(Protocol):
46 """Operations a parametric active-set path tracer needs from its Hessian.
48 The turning-point loop touches the symmetric positive (semi-)definite matrix
49 ``H`` of the quadratic objective through a small number of operations: full
50 matrix-vector products, solves against the *free* (active) principal block,
51 and cross-products between free and blocked coordinates. In the Critical Line
52 Algorithm ``H`` is the covariance ``Sigma``; in a LASSO / LARS path it is the
53 Gram matrix ``X.T @ X``. Any object implementing this protocol can serve as
54 that backend.
56 The reference implementation is ``DenseCovariance``, which wraps an explicit
57 matrix; structured backends (e.g. diagonal-plus-low-rank via the Woodbury
58 identity) implement the same contract without ever materialising an
59 ``n x n`` matrix.
61 ``CovarianceOperator`` is kept as a backward-compatible alias of this
62 protocol (see ``CovarianceOperator`` below).
63 """
65 @property
66 def n(self) -> int:
67 """Number of variables (the dimension of the quadratic form)."""
68 ... # pragma: no cover
70 def matvec(self, x: NDArray[np.float64]) -> NDArray[np.float64]:
71 """Compute the matrix-vector product ``Sigma @ x``.
73 Args:
74 x: Vector of shape ``(n,)`` or matrix of shape ``(n, r)``.
76 Returns:
77 ``Sigma @ x`` with the same trailing shape as ``x``.
78 """
79 ... # pragma: no cover
81 def solve_free(self, free: NDArray[np.bool_], rhs: NDArray[np.float64]) -> NDArray[np.float64]:
82 """Solve the free-block system ``Sigma[free][:, free] @ y = rhs``.
84 Args:
85 free: Boolean mask of shape ``(n,)`` selecting the free assets.
86 rhs: Right-hand side of shape ``(n_free,)`` or ``(n_free, r)``.
88 Returns:
89 The solution ``y`` with the same shape as ``rhs``.
90 """
91 ... # pragma: no cover
93 def cross(self, free: NDArray[np.bool_], x: NDArray[np.float64]) -> NDArray[np.float64]:
94 """Compute the free-to-blocked cross product ``Sigma[free][:, ~free] @ x[~free]``.
96 Args:
97 free: Boolean mask of shape ``(n,)`` selecting the free assets.
98 x: Full-length vector of shape ``(n,)``; only the blocked entries
99 ``x[~free]`` enter the product.
101 Returns:
102 Vector of shape ``(n_free,)``.
103 """
104 ... # pragma: no cover
106 def rcond_free(self, free: NDArray[np.bool_]) -> float:
107 """Reciprocal condition number of the free block ``Sigma[free][:, free]``.
109 Returns a value in ``[0, 1]`` (``1`` well conditioned, ``0`` singular).
110 The CLA uses this to detect a rank-deficient free block, whose solve is
111 unreliable, directly and portably, rather than inferring it from the
112 (backend-dependent) magnitude of the resulting box violation.
114 Args:
115 free: Boolean mask of shape ``(n,)`` selecting the free assets.
117 Returns:
118 The reciprocal 2-norm condition number of the free block.
119 """
120 ... # pragma: no cover
123def bordered_solve(
124 quad: QuadraticForm,
125 free: NDArray[np.bool_],
126 c_free: NDArray[np.float64],
127 rhs_const: NDArray[np.float64],
128 rhs_slope: NDArray[np.float64],
129 d_const: NDArray[np.float64],
130 d_slope: NDArray[np.float64],
131) -> tuple[
132 NDArray[np.float64],
133 NDArray[np.float64],
134 NDArray[np.float64],
135 NDArray[np.float64],
136]:
137 """Solve the bordered KKT system ``[[H_FF, C_F^T], [C_F, 0]] @ [x; nu] = [rhs; d]``.
139 The shared segment solve of ``cvxcla.cla`` and ``cvxcla.lasso``. ``H_FF`` is the
140 free block of the quadratic form, ``C_F`` the constraint rows on the free columns.
141 The constant and ``lambda``-slope right-hand sides are solved together so ``H_FF``
142 is factorised once; ``mc == 0`` (no constraint rows) is the plain free-block solve.
143 ``H_FF`` is reached only through ``solve_free``, so structured backends never
144 materialise an ``n x n`` matrix. Returns ``(x_const, x_slope, nu_const, nu_slope)``:
145 the free solution and constraint multipliers per system (multipliers empty if
146 ``mc == 0``). Callers assemble the right-hand sides and interpret the outputs.
147 """
148 mc = c_free.shape[0]
149 if mc == 0:
150 sol = quad.solve_free(free, np.column_stack([rhs_const, rhs_slope]))
151 empty: NDArray[np.float64] = np.zeros(0)
152 return sol[:, 0], sol[:, 1], empty, empty
154 # One multi-RHS solve against H_FF covers the constraint columns C_F^T and both
155 # the constant and slope systems, so H_FF is factorised a single time.
156 solved = quad.solve_free(free, np.column_stack([c_free.T, rhs_const, rhs_slope]))
157 y = solved[:, :mc] # H_FF^{-1} C_F^T
158 z_const = solved[:, mc]
159 z_slope = solved[:, mc + 1]
161 # Schur complement C_F H_FF^{-1} C_F^T and the stacked multipliers.
162 schur = c_free @ y
163 nu = cast(
164 "NDArray[np.float64]",
165 np.linalg.solve(schur, np.column_stack([c_free @ z_const - d_const, c_free @ z_slope - d_slope])),
166 )
167 nu_const, nu_slope = nu[:, 0], nu[:, 1]
169 # Back-substitute the free solution.
170 return z_const - y @ nu_const, z_slope - y @ nu_slope, nu_const, nu_slope
173# Backward-compatible alias. The protocol was introduced as ``CovarianceOperator``
174# for the portfolio (covariance) setting; it is the Hessian contract for any
175# parametric active-set path, so the canonical name is now ``QuadraticForm``.
176# Existing imports of ``CovarianceOperator`` keep working unchanged.
177CovarianceOperator = QuadraticForm