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

1"""Feasibility projection for Critical Line Algorithm turning points. 

2 

3Accumulated floating-point round-off over the many turning points of a long 

4trace can place a candidate a hair outside its box even at a well-conditioned 

5vertex: the covariance there has near-flat directions (its small eigenvalues) 

6and the round-off lies in exactly those directions, so the candidate is optimal 

7to solver precision but not exactly feasible. These pure helpers project such a 

8candidate back onto the feasible region ``{w : lower <= w <= upper, C w = d}``. 

9 

10They are a strict no-op for the well-posed, already-feasible turning points that 

11are the common case, so they never perturb the exact frontier. ``CLA._emit`` 

12calls :func:`project_feasible` at every turning point. 

13""" 

14 

15from __future__ import annotations 

16 

17import numpy as np 

18from cvx.linalg import AffineProjection 

19from numpy.typing import NDArray 

20 

21 

22def project_feasible( 

23 weights: NDArray[np.float64], 

24 lower: NDArray[np.float64], 

25 upper: NDArray[np.float64], 

26 a: NDArray[np.float64], 

27 b: NDArray[np.float64], 

28 g: NDArray[np.float64], 

29 h: NDArray[np.float64], 

30 active_ineq: NDArray[np.bool_], 

31) -> NDArray[np.float64]: 

32 """Project ``weights`` onto the feasible region, clearing round-off. 

33 

34 Well-posed turning points are already strictly feasible and are returned 

35 unchanged so the projection never perturbs the exact frontier. Otherwise the 

36 candidate is dispatched to the closed-form capped-simplex projection for the 

37 canonical all-ones budget with no active inequality row (see 

38 :func:`project_capped_simplex`), and to the general alternating projection 

39 otherwise (see :func:`project_alternating`). 

40 

41 Args: 

42 weights: The candidate weight vector to project. 

43 lower: Per-asset lower bounds. 

44 upper: Per-asset upper bounds. 

45 a: Equality-constraint matrix ``A`` of ``A w = b``. 

46 b: Equality-constraint right-hand side ``b``. 

47 g: Inequality-constraint matrix ``G`` of ``G w <= h`` (``(p, n)``). 

48 h: Inequality-constraint right-hand side ``h`` (length ``p``). 

49 active_ineq: Boolean mask (length ``p``) of the active inequality rows. 

50 

51 Returns: 

52 The projected weight vector, feasible to the box and the constraints. 

53 """ 

54 # Well-posed turning points are already strictly feasible: return them 

55 # unchanged so the projection never perturbs the exact frontier. 

56 if np.all(weights >= lower) and np.all(weights <= upper): 

57 return weights 

58 

59 if _is_capped_simplex(active_ineq, a): 

60 return project_capped_simplex(weights, lower, upper, float(b[0])) 

61 

62 c = np.vstack([a, g[active_ineq]]) 

63 d = np.concatenate([b, h[active_ineq]]) 

64 return project_alternating(weights, lower, upper, c, d) 

65 

66 

67def _is_capped_simplex(active_ineq: NDArray[np.bool_], a: NDArray[np.float64]) -> bool: 

68 """Whether the projection reduces to the closed-form capped simplex. 

69 

70 True when no inequality row is active and the equality is the canonical 

71 all-ones budget (a single row of ones), the case handled by 

72 :func:`project_capped_simplex`. 

73 

74 Args: 

75 active_ineq: Boolean mask of the active inequality rows. 

76 a: Equality-constraint matrix ``A``. 

77 

78 Returns: 

79 ``True`` if the closed-form capped-simplex projection applies. 

80 """ 

81 return not active_ineq.any() and a.shape[0] == 1 and bool(np.allclose(a, 1.0)) 

82 

83 

84def project_capped_simplex( 

85 weights: NDArray[np.float64], 

86 lower: NDArray[np.float64], 

87 upper: NDArray[np.float64], 

88 budget: float, 

89) -> NDArray[np.float64]: 

90 """Euclidean projection onto the capped simplex for the all-ones budget. 

91 

92 Computed by water-filling a single shift ``theta`` so that 

93 ``sum(clip(w - theta, lower, upper)) = budget``. A plain clip-then-rescale is 

94 deliberately *not* used: rescaling the clipped weights to restore the budget 

95 can push capped weights back over their bound when many assets are capped at 

96 once (heavy ties under a tight cap), re-introducing the very infeasibility 

97 the projection is meant to clear. 

98 

99 Args: 

100 weights: The candidate weight vector, known to violate its box. 

101 lower: Per-asset lower bounds. 

102 upper: Per-asset upper bounds. 

103 budget: The all-ones budget right-hand side ``sum(w)``. 

104 

105 Returns: 

106 The projected weight vector, on the budget and inside the box. 

107 """ 

108 # sum(clip(w - theta, lower, upper)) is non-increasing in theta; bisect 

109 # for the theta that hits the budget. The bracket clips to all-upper 

110 # (sum at least the budget) at theta_lo and all-lower (at most the 

111 # budget) at theta_hi, so a root is guaranteed for any feasible problem. 

112 theta_lo = float((weights - upper).min()) - 1.0 

113 theta_hi = float((weights - lower).max()) + 1.0 

114 for _ in range(100): 

115 theta = 0.5 * (theta_lo + theta_hi) 

116 if float(np.clip(weights - theta, lower, upper).sum()) > budget: 

117 theta_lo = theta 

118 else: 

119 theta_hi = theta 

120 return np.clip(weights - 0.5 * (theta_lo + theta_hi), lower, upper) 

121 

122 

123def project_alternating( 

124 weights: NDArray[np.float64], 

125 lower: NDArray[np.float64], 

126 upper: NDArray[np.float64], 

127 c: NDArray[np.float64], 

128 d: NDArray[np.float64], 

129) -> NDArray[np.float64]: 

130 """Alternating projection onto the box and the affine set ``{C w = d}``. 

131 

132 ``C`` stacks the equality rows ``A`` and the active inequality rows ``G_S`` 

133 (held at equality ``g_i w = h_i``). The candidate already satisfies 

134 ``C w = d`` (the reduced KKT solve enforces it) and the inactive inequality 

135 rows keep a margin, so a few iterations alternating a box clip with the 

136 affine projection converge to a point feasible to the box, the equalities, 

137 and every inequality. 

138 

139 Args: 

140 weights: The candidate weight vector, known to violate its box. 

141 lower: Per-asset lower bounds. 

142 upper: Per-asset upper bounds. 

143 c: Stacked equality/active-inequality matrix. 

144 d: Stacked equality/active-inequality right-hand side. 

145 

146 Returns: 

147 The projected weight vector, feasible to the box and the constraints. 

148 """ 

149 affine = AffineProjection(c, d) 

150 projected = weights 

151 for _ in range(100): 

152 projected = np.clip(projected, lower, upper) 

153 projected = affine.project(projected) 

154 return np.clip(projected, lower, upper)