Coverage for src/cvxcla/operators/_core.py: 100%

1"""Operator protocol alias and the parametric-path helpers built on cvx-linalg. 

2 

3The parametric active-set path tracer reaches its Hessian through the cvx-linalg 

4symmetric-operator protocol (``matvec`` / ``block_matvec`` / ``solve_free`` / 

5``rcond_free``). :data:`QuadraticForm` is that contract; :data:`CovarianceOperator` 

6is a backward-compatible alias for the portfolio (covariance) setting. The 

7concrete backends live in :mod:`cvx.linalg`; :mod:`cvxcla.operators.builders` 

8assembles them from CLA / LASSO inputs. 

9 

10The generic linear algebra -- the bordered KKT (Schur complement) solve and the 

11affine projection -- also lives in :mod:`cvx.linalg`. What remains here is the 

12thin homotopy-specific glue: adapting the loop's boolean masks to the operators' 

13integer-index API, and packing the constant / ``lambda``-slope pair of a 

14parametric segment into the shared multi-RHS solve. 

15""" 

16 

17from __future__ import annotations 

18 

19import numpy as np 

20from cvx.linalg import SymmetricOperator 

21from cvx.linalg import bordered_solve as _bordered_solve 

22from numpy.typing import NDArray 

23 

24# The Hessian contract for a parametric active-set path. In the CLA it is the 

25# covariance ``Sigma``; in a LASSO / LARS path it is the Gram matrix ``X.T @ X``. 

26QuadraticForm = SymmetricOperator 

27CovarianceOperator = SymmetricOperator 

28 

29# The singularity threshold for ``rcond_free``. A genuinely rank-deficient free 

30# block has a reciprocal condition number at round-off level (~1e-16); a 

31# well-posed or merely near-degenerate block sits many orders above it (>= ~1e-4 

32# across the degeneracy sweep in experiments/degeneracy_boundary.py). The 1e-12 

33# cut sits in the wide gap between the two and is the conventional 

34# numerical-singularity scale. 

35_RCOND_FLOOR = 1e-12 # pragma: no mutate 

36 

37 

38def cross(operator: SymmetricOperator, free: NDArray[np.bool_], x: NDArray[np.float64]) -> NDArray[np.float64]: 

39 """Free-to-blocked cross product ``H[free][:, ~free] @ x[~free]`` from a boolean mask. 

40 

41 Adapts the boolean-mask indexing the path tracers use to the integer-index 

42 :meth:`~cvx.linalg.SymmetricOperator.block_matvec` of a symmetric operator. 

43 

44 Args: 

45 operator: The symmetric operator (Hessian) backend. 

46 free: Boolean mask of shape ``(n,)`` selecting the free coordinates. 

47 x: Full-length vector of shape ``(n,)``; only ``x[~free]`` enters the product. 

48 

49 Returns: 

50 Vector of shape ``(n_free,)``. 

51 """ 

52 result: NDArray[np.float64] = operator.block_matvec(np.flatnonzero(free), np.flatnonzero(~free), x[~free]) 

53 return result 

54 

55 

56def bordered_solve( 

57 quad: SymmetricOperator, 

58 free: NDArray[np.bool_], 

59 c_free: NDArray[np.float64], 

60 rhs_const: NDArray[np.float64], 

61 rhs_slope: NDArray[np.float64], 

62 d_const: NDArray[np.float64], 

63 d_slope: NDArray[np.float64], 

64) -> tuple[ 

65 NDArray[np.float64], 

66 NDArray[np.float64], 

67 NDArray[np.float64], 

68 NDArray[np.float64], 

69]: 

70 """Solve the bordered KKT system for a parametric segment's constant and slope parts. 

71 

72 A thin adapter over :func:`cvx.linalg.bordered_solve`: it converts the loop's 

73 boolean *free* mask to integer indices and packs the constant and ``lambda``-slope 

74 right-hand sides as the two columns of one multi-RHS solve, so ``H_FF`` and the 

75 Schur complement are factorised once. Returns 

76 ``(x_const, x_slope, nu_const, nu_slope)`` (multipliers empty when there are no 

77 constraint rows). 

78 """ 

79 x, nu = _bordered_solve( 

80 quad, 

81 np.flatnonzero(free), 

82 c_free, 

83 np.column_stack([rhs_const, rhs_slope]), 

84 np.column_stack([d_const, d_slope]), 

85 ) 

86 return x[:, 0], x[:, 1], nu[:, 0], nu[:, 1]