Coverage for src/fast_minimum_variance/_base.py: 99%

1"""Common base for portfolio-optimisation problem classes.""" 

2 

3from abc import ABC, abstractmethod 

4from collections.abc import Callable 

5from dataclasses import dataclass 

6from typing import Any 

7 

8import cvxpy as cp 

9import numpy as np 

10from cvx.linalg import cholesky 

11 

12 

13@dataclass(frozen=True) 

14class _BaseProblem(ABC): 

15 """Shared fields, utilities, and solver templates for portfolio problems. 

16 

17 Subclasses must implement the four abstract hooks: 

18 

19 * ``_constraint_active_set(solve_fn)`` — outer constraint-handling loop 

20 * ``_kkt_step(mask) -> (w, iters)`` — one direct-KKT inner step 

21 * ``_cg_step(mask) -> (w, iters)`` — one CG inner step 

22 * ``_cvxpy_constraints(w, cp) -> list`` — CVXPY constraint list 

23 

24 All ``solve_*`` methods are implemented here as template methods that 

25 call ``_constraint_active_set`` with the appropriate ``_XXX_step`` 

26 method, then optionally clip-and-renormalize. 

27 """ 

28 

29 X: np.ndarray 

30 target: np.ndarray | None = None 

31 alpha: float = 0.0 

32 rho: float = 0.0 

33 mu: np.ndarray | None = None 

34 target_lr: tuple[float, np.ndarray, np.ndarray] | None = None # (bar_lam, U_k, delta_k) — low-rank + identity 

35 pcg_lr: tuple[float, np.ndarray, np.ndarray] | None = None # (bar_lam, U_k, delta_k) — RMT preconditioner (§5.3) 

36 

37 def __post_init__(self) -> None: 

38 """Validate target/target_lr shapes when supplied.""" 

39 n = self.n 

40 if self.target is not None and self.target.shape != (n, n): 

41 raise ValueError(f"target must be a square {n} x {n} matrix, got {self.target.shape}") # noqa: TRY003 

42 if self.target_lr is not None: 

43 _bar_lam, U_k, delta_k = self.target_lr # noqa: N806 

44 if U_k.shape[0] != n or U_k.shape[1] != delta_k.shape[0]: 

45 raise ValueError( # noqa: TRY003 

46 f"target_lr: U_k must be ({n}, k) and delta_k (k,), got {U_k.shape}, {delta_k.shape}" 

47 ) 

48 

49 # ------------------------------------------------------------------ 

50 # Shared utilities 

51 # ------------------------------------------------------------------ 

52 

53 @property 

54 def t(self) -> int: 

55 """Return the number of rows in X.""" 

56 return int(self.X.shape[0]) 

57 

58 @property 

59 def n(self) -> int: 

60 """Number of assets (columns of X).""" 

61 return int(self.X.shape[1]) 

62 

63 @staticmethod 

64 def _clip_and_renormalize(w: np.ndarray) -> np.ndarray: 

65 """Clip weights to ``[0, ∞)`` and renormalize to sum to 1.""" 

66 w = np.maximum(w, 0) 

67 w /= w.sum() 

68 return w 

69 

70 # ------------------------------------------------------------------ 

71 # Abstract hooks (raise NotImplementedError — subclasses must override) 

72 # ------------------------------------------------------------------ 

73 @abstractmethod 

74 def _constraint_active_set( 

75 self, 

76 solve_fn: Callable[[np.ndarray], tuple[np.ndarray, int]], 

77 tol: float = 1e-6, 

78 max_iter: int = 10_000, 

79 ) -> tuple[np.ndarray, int, int]: # pragma: no cover 

80 """Run the outer constraint-handling loop, calling ``solve_fn`` each iteration.""" 

81 raise NotImplementedError 

82 

83 @abstractmethod 

84 def _kkt_step(self, active: np.ndarray) -> tuple[np.ndarray, int]: # pragma: no cover 

85 """Solve one inner direct-KKT step; return ``(w, iters)``.""" 

86 raise NotImplementedError 

87 

88 @abstractmethod 

89 def _cvxpy_constraints(self, w: cp.Variable, cp: object) -> list[Any]: # pragma: no cover 

90 """Return the list of CVXPY constraints for ``solve_cvxpy``.""" 

91 raise NotImplementedError 

92 

93 @abstractmethod 

94 def _cg_step(self, active: np.ndarray) -> tuple[np.ndarray, int]: 

95 """Solve one inner CG step; return ``(w, iters)``.""" 

96 raise NotImplementedError # pragma: no cover 

97 

98 def _pcg_step(self, active: np.ndarray, x0: np.ndarray | None = None) -> tuple[np.ndarray, int]: # pragma: no cover 

99 """Solve one inner PCG step with RMT preconditioner; return ``(w, iters)``. 

100 

101 Subclasses that support PCG (e.g. ``_MinVarProblem``) override this. 

102 The base implementation raises so callers get a clear error if PCG is 

103 invoked on a problem type that has not implemented it. 

104 """ 

105 raise NotImplementedError 

106 

107 # ------------------------------------------------------------------ 

108 # Template solvers 

109 # ------------------------------------------------------------------ 

110 

111 def solve_kkt(self, *, project: bool = True) -> tuple[np.ndarray, int]: 

112 """Solve via the direct KKT system. 

113 

114 Args: 

115 project: Clip weights to ``[0, ∞)`` and renormalize to sum to 1 

116 after solving. Set to ``False`` for custom constraints. 

117 

118 Returns: 

119 ``(w, n_iters)`` — weight vector of shape ``(N,)`` and number of 

120 outer iterations taken. 

121 

122 Examples: 

123 >>> import numpy as np 

124 >>> from fast_minimum_variance import Problem 

125 >>> X = np.random.default_rng(0).standard_normal((100, 5)) 

126 >>> w, iters = Problem(X).solve_kkt() 

127 >>> float(round(w.sum(), 10)) 

128 1.0 

129 >>> bool((w >= 0).all()) 

130 True 

131 """ 

132 w, outer, _inner = self._constraint_active_set(self._kkt_step) 

133 if project: 

134 w = self._clip_and_renormalize(w) 

135 return w, outer 

136 

137 def solve_cvxpy(self, *, project: bool = True) -> tuple[np.ndarray, int]: 

138 """Solve via CVXPY with the Clarabel backend (reference solver). 

139 

140 Requires the ``convex`` extra:: 

141 

142 pip install fast-minimum-variance[convex] 

143 

144 Args: 

145 project: Clip and renormalize after solving (see ``solve_kkt``). 

146 

147 Returns: 

148 ``(w, n_iters)`` — weight vector of shape ``(N,)`` and solver 

149 iteration count. 

150 

151 Examples: 

152 >>> import numpy as np 

153 >>> from fast_minimum_variance import Problem 

154 >>> X = np.random.default_rng(0).standard_normal((100, 5)) 

155 >>> w, iters = Problem(X).solve_cvxpy() 

156 >>> float(round(w.sum(), 6)) 

157 1.0 

158 >>> bool((w >= -1e-6).all()) 

159 True 

160 """ 

161 w = cp.Variable(self.n) 

162 if self.target is not None: 

163 # target is the penalty matrix M; decompose as M = chol chol^T so ||chol^T w||^2 = w^T M w 

164 chol = cholesky(self.target) 

165 objective = (1.0 - self.alpha) * cp.sum_squares(self.X @ w) / self.t + self.alpha * cp.sum_squares( 

166 chol.T @ w 

167 ) 

168 else: 

169 objective = cp.sum_squares(self.X @ w) / self.t 

170 if self.rho != 0.0 and self.mu is not None: 

171 objective = objective - self.rho * (self.mu @ w) 

172 

173 problem = cp.Problem(cp.Minimize(objective), self._cvxpy_constraints(w, cp)) 

174 problem.solve(solver=cp.CLARABEL) 

175 

176 result = w.value 

177 if result is None: 

178 raise RuntimeError("CVXPY solver failed to find a solution") # noqa: TRY003 

179 if project: 

180 result = self._clip_and_renormalize(result) 

181 return result, int(problem.solver_stats.num_iters or 0) 

182 

183 def solve_cg(self, *, project: bool = True) -> tuple[np.ndarray, int, int]: 

184 """Solve via matrix-free conjugate gradients. 

185 

186 Args: 

187 project: Clip weights to ``[0, ∞)`` and renormalize to sum to 1 

188 after solving. Set to ``False`` for custom constraints. 

189 

190 Returns: 

191 ``(w, outer_steps, inner_iters)`` — weight vector, number of outer 

192 active-set steps, and total CG iterations summed across all steps. 

193 

194 Examples: 

195 >>> import numpy as np 

196 >>> from fast_minimum_variance import Problem 

197 >>> X = np.random.default_rng(0).standard_normal((100, 5)) 

198 >>> w, outer, inner = Problem(X).solve_cg() 

199 >>> float(round(w.sum(), 10)) 

200 1.0 

201 >>> bool((w >= 0).all()) 

202 True 

203 """ 

204 w, outer, inner = self._constraint_active_set(self._cg_step) 

205 if project: 

206 w = self._clip_and_renormalize(w) 

207 return w, outer, inner 

208 

209 def solve_pcg(self, *, project: bool = True) -> tuple[np.ndarray, int, int]: 

210 """Solve via matrix-free PCG with RMT preconditioner (Section 5.3). 

211 

212 Solves ``Sigma_LW_oracle x = 1`` using ``T0^RMT`` as preconditioner. 

213 Requires ``pcg_lr = (bar_lam, U_k, delta_k)`` from RMT preprocessing. 

214 The preconditioner is applied via the Woodbury identity at O(nk) per step; 

215 the system matvec costs O(nT). Returns the oracle-LW minimum-variance 

216 portfolio — not the RMT portfolio — in O(sqrt(1/alpha_oracle)) iterations. 

217 

218 Returns: 

219 ``(w, outer_steps, inner_iters)`` 

220 

221 Examples: 

222 >>> import numpy as np 

223 >>> from fast_minimum_variance import Problem 

224 >>> rng = np.random.default_rng(0) 

225 >>> X = rng.standard_normal((100, 5)) 

226 >>> bar_lam = float(np.trace(X.T @ X / 100) / 5) 

227 >>> U_k = np.eye(5, 2) 

228 >>> delta_k = np.array([0.1, 0.05]) 

229 >>> w, outer, inner = Problem(X, alpha=0.1, pcg_lr=(bar_lam, U_k, delta_k)).solve_pcg() 

230 >>> float(round(w.sum(), 10)) 

231 1.0 

232 >>> bool((w >= 0).all()) 

233 True 

234 """ 

235 if self.pcg_lr is None: 

236 raise ValueError("pcg_lr must be set; pass pcg_lr=(bar_lam, U_k, delta_k)") # noqa: TRY003 

237 w, outer, inner = self._constraint_active_set(self._pcg_step) 

238 if project: 

239 w = self._clip_and_renormalize(w) 

240 return w, outer, inner