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

1"""Builders that assemble cvx-linalg symmetric operators from CLA / LASSO inputs. 

2 

3These replace the former operator *classes*. The numerics -- matrix-vector and 

4block products, free-block solves (Cholesky, Woodbury, maintained inverse), and 

5the reciprocal-condition check -- now live in :mod:`cvx.linalg`. Here we only 

6build the right operator from the domain inputs: a covariance matrix, a 

7returns / data matrix, or a factor model. The backward-compatible names 

8``DenseCovariance`` / ``IncrementalDenseCovariance`` / ``GramCovariance`` / 

9``FactorCovariance`` are kept as aliases of these builders. 

10""" 

11 

12from __future__ import annotations 

13 

14import numpy as np 

15from cvx.linalg import DenseOperator, FactorOperator, GramOperator, IncrementalDenseOperator 

16from numpy.typing import NDArray 

17 

18 

19def _symmetric_matrix(matrix: NDArray[np.float64]) -> NDArray[np.float64]: 

20 """Return *matrix* as a float array after checking it is square and symmetric. 

21 

22 Raises: 

23 ValueError: If *matrix* is not square or not symmetric to tolerance. 

24 """ 

25 matrix = np.asarray(matrix, dtype=np.float64) 

26 if matrix.ndim != 2 or matrix.shape[0] != matrix.shape[1]: 

27 msg = f"Covariance must be a square matrix, got shape {matrix.shape}" 

28 raise ValueError(msg) 

29 if not np.allclose(matrix, matrix.T): 

30 msg = "Covariance must be symmetric" 

31 raise ValueError(msg) 

32 return matrix 

33 

34 

35def dense_covariance(matrix: NDArray[np.float64]) -> DenseOperator: 

36 """Build a :class:`~cvx.linalg.DenseOperator` from an explicit symmetric covariance. 

37 

38 Args: 

39 matrix: A symmetric ``(n, n)`` covariance matrix. 

40 

41 Returns: 

42 A dense symmetric operator wrapping *matrix*. 

43 

44 Examples: 

45 >>> import numpy as np 

46 >>> from cvxcla import DenseCovariance 

47 >>> sigma = np.array([[4.0, 1.0], [1.0, 9.0]]) 

48 >>> op = DenseCovariance(sigma) 

49 >>> op.n 

50 2 

51 >>> op.matvec(np.array([1.0, 0.0])) 

52 array([4., 1.]) 

53 

54 A non-square or non-symmetric input is rejected rather than silently 

55 symmetrised, so a transposed or half-filled matrix fails at construction 

56 instead of producing a wrong frontier: 

57 

58 >>> DenseCovariance(np.array([[1.0, 2.0], [0.0, 1.0]])) 

59 Traceback (most recent call last): 

60 ... 

61 ValueError: Covariance must be symmetric 

62 """ 

63 return DenseOperator(_symmetric_matrix(matrix)) 

64 

65 

66def incremental_dense_covariance(matrix: NDArray[np.float64]) -> IncrementalDenseOperator: 

67 """Build an :class:`~cvx.linalg.IncrementalDenseOperator` (maintained free-block inverse). 

68 

69 A drop-in alternative to :func:`dense_covariance` for a loop that changes its 

70 free set one index at a time; see the operator's own caveats on numerics. 

71 

72 Args: 

73 matrix: A symmetric ``(n, n)`` covariance matrix. 

74 

75 Returns: 

76 A dense symmetric operator that maintains the free-block inverse. 

77 

78 Examples: 

79 >>> import numpy as np 

80 >>> from cvxcla import CLA, DenseCovariance, IncrementalDenseCovariance 

81 >>> sigma = np.array([[4.0, 1.0], [1.0, 9.0]]) 

82 >>> x = np.array([1.0, 2.0]) 

83 >>> incremental = IncrementalDenseCovariance(sigma) 

84 >>> bool(np.allclose(incremental.matvec(x), DenseCovariance(sigma).matvec(x))) 

85 True 

86 

87 It is a drop-in swap: the two agree on the same problem, and only the 

88 cost profile of the free-block solves differs. 

89 

90 >>> cla = CLA.problem(np.array([0.1, 0.2]), incremental).long_only().budget().trace() 

91 >>> len(cla) > 0 

92 True 

93 """ 

94 return IncrementalDenseOperator(_symmetric_matrix(matrix)) 

95 

96 

97def gram_covariance(returns: NDArray[np.float64], ridge: float = 0.0) -> GramOperator: 

98 """Build a :class:`~cvx.linalg.GramOperator` for the sample covariance of *returns*. 

99 

100 The sample covariance ``X_c.T X_c / (T - 1)`` (``X_c`` the column-centered 

101 ``(T, n)`` data) is realised by absorbing the scale into the factor: 

102 ``M = sqrt(1 / (T - 1)) * X_c``, so the operator represents 

103 ``M.T M + ridge * I`` and never forms the ``n x n`` covariance. 

104 

105 Args: 

106 returns: The ``(T, n)`` matrix of observations (``T >= 2``). 

107 ridge: A non-negative diagonal loading added to the covariance. 

108 

109 Returns: 

110 A matrix-free Gram operator for the (ridged) sample covariance. 

111 

112 Raises: 

113 ValueError: If *returns* is not a ``(T, n)`` matrix with ``T >= 2``. 

114 

115 Examples: 

116 >>> import numpy as np 

117 >>> from cvxcla import GramCovariance 

118 >>> returns = np.array([[1.0, 2.0], [3.0, 5.0], [5.0, 4.0]]) 

119 >>> op = GramCovariance(returns) 

120 

121 The operator represents exactly the sample covariance ``np.cov`` would 

122 form -- without ever building the ``(n, n)`` matrix: 

123 

124 >>> sample = np.cov(returns, rowvar=False) 

125 >>> bool(np.allclose(op.matvec(np.array([1.0, 2.0])), sample @ np.array([1.0, 2.0]))) 

126 True 

127 >>> bool(np.allclose(op.diag, sample.diagonal())) 

128 True 

129 

130 ``ridge`` adds a diagonal loading, which is how a short sample (fewer 

131 observations than assets) is made positive definite: 

132 

133 >>> float(np.round(GramCovariance(returns, ridge=0.5).diag[0] - op.diag[0], 12)) 

134 0.5 

135 

136 A single observation cannot define a covariance and is refused: 

137 

138 >>> GramCovariance(np.array([[1.0, 2.0]])) 

139 Traceback (most recent call last): 

140 ... 

141 ValueError: returns must be a (T, n) matrix with T >= 2 observations, got shape (1, 2) 

142 """ 

143 returns = np.asarray(returns, dtype=np.float64) 

144 if returns.ndim != 2 or returns.shape[0] < 2: 

145 msg = f"returns must be a (T, n) matrix with T >= 2 observations, got shape {returns.shape}" 

146 raise ValueError(msg) 

147 t = returns.shape[0] 

148 centered = returns - returns.mean(axis=0, keepdims=True) 

149 factor = centered / np.sqrt(t - 1.0) 

150 return GramOperator(factor, ridge=ridge) 

151 

152 

153def factor_covariance( 

154 d: NDArray[np.float64], 

155 u: NDArray[np.float64], 

156 delta: NDArray[np.float64], 

157) -> FactorOperator: 

158 """Build a :class:`~cvx.linalg.FactorOperator` for ``Sigma = diag(d) + U Delta U.T``. 

159 

160 Args: 

161 d: Positive idiosyncratic variances of shape ``(n,)``. 

162 u: Factor loadings of shape ``(n, k)``. 

163 delta: Factor covariance, either ``(k,)`` eigenvalues (a diagonal ``Delta``) 

164 or a symmetric positive-definite ``(k, k)`` matrix. 

165 

166 Returns: 

167 A diagonal-plus-low-rank operator with Woodbury free-block solves. 

168 

169 Raises: 

170 ValueError: If *delta* is neither a ``(k,)`` vector nor a ``(k, k)`` matrix. 

171 

172 Examples: 

173 >>> import numpy as np 

174 >>> from cvxcla import FactorCovariance 

175 >>> d = np.array([0.5, 0.5, 0.5]) 

176 >>> u = np.array([[1.0], [1.0], [-1.0]]) 

177 >>> delta = np.array([2.0]) 

178 >>> op = FactorCovariance(d=d, u=u, delta=delta) 

179 >>> op.n, op.k 

180 (3, 1) 

181 

182 The operator matches its dense materialisation while storing only 

183 ``d`` and ``u`` -- the point of the diagonal-plus-low-rank form: 

184 

185 >>> dense = np.diag(d) + u @ np.diag(delta) @ u.T 

186 >>> x = np.array([1.0, 0.0, 0.0]) 

187 >>> bool(np.allclose(op.matvec(x), dense @ x)) 

188 True 

189 

190 ``delta`` may equally be a full ``(k, k)`` factor covariance: 

191 

192 >>> bool(np.allclose(FactorCovariance(d=d, u=u, delta=np.diag(delta)).matvec(x), dense @ x)) 

193 True 

194 

195 Anything else is refused: 

196 

197 >>> FactorCovariance(d=d, u=u, delta=np.zeros((1, 1, 1))) 

198 Traceback (most recent call last): 

199 ... 

200 ValueError: delta must be a (k,) vector or (k, k) matrix, got ndim 3 

201 """ 

202 d = np.asarray(d, dtype=np.float64) 

203 u = np.asarray(u, dtype=np.float64) 

204 delta = np.asarray(delta, dtype=np.float64) 

205 if delta.ndim == 1: 

206 inner = np.diag(delta) 

207 elif delta.ndim == 2: 

208 inner = delta 

209 else: 

210 msg = f"delta must be a (k,) vector or (k, k) matrix, got ndim {delta.ndim}" 

211 raise ValueError(msg) 

212 return FactorOperator(d, u, inner) 

213 

214 

215# Backward-compatible names: the operator *classes* are gone, but the familiar 

216# constructor-style names remain as builders returning cvx-linalg operators. 

217DenseCovariance = dense_covariance 

218IncrementalDenseCovariance = incremental_dense_covariance 

219GramCovariance = gram_covariance 

220FactorCovariance = factor_covariance