Coverage for src / cvx / risk / rand / rand_cov.py: 100%

1"""Random covariance matrix generation utilities. 

2 

3This module provides functions for generating random positive semi-definite 

4covariance matrices. These are useful for testing and simulation purposes. 

5 

6Example: 

7 Generate a random covariance matrix: 

8 

9 >>> import numpy as np 

10 >>> from cvx.risk.rand import rand_cov 

11 >>> # Generate a 5x5 random covariance matrix 

12 >>> cov = rand_cov(5, seed=42) 

13 >>> cov.shape 

14 (5, 5) 

15 >>> # Verify it's symmetric 

16 >>> bool(np.allclose(cov, cov.T)) 

17 True 

18 >>> # Verify it's positive semi-definite 

19 >>> bool(np.all(np.linalg.eigvals(cov) >= -1e-10)) 

20 True 

21 

22""" 

23 

24# Copyright 2023 Stanford University Convex Optimization Group 

25# 

26# Licensed under the Apache License, Version 2.0 (the "License"); 

27# you may not use this file except in compliance with the License. 

28# You may obtain a copy of the License at 

29# 

30# http://www.apache.org/licenses/LICENSE-2.0 

31# 

32# Unless required by applicable law or agreed to in writing, software 

33# distributed under the License is distributed on an "AS IS" BASIS, 

34# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 

35# See the License for the specific language governing permissions and 

36# limitations under the License. 

37from __future__ import annotations 

38 

39import numpy as np 

40 

41 

42def rand_cov(n: int, seed: int | None = None) -> np.ndarray: 

43 """Construct a random positive semi-definite covariance matrix of size n x n. 

44 

45 The matrix is constructed as A^T @ A where A is a random n x n matrix with 

46 elements drawn from a standard normal distribution. This ensures the result 

47 is symmetric and positive semi-definite. 

48 

49 Args: 

50 n: Size of the covariance matrix (n x n). 

51 seed: Random seed for reproducibility. If None, uses the current 

52 random state. 

53 

54 Returns: 

55 A random positive semi-definite n x n covariance matrix. 

56 

57 Example: 

58 Generate a reproducible random covariance matrix: 

59 

60 >>> import numpy as np 

61 >>> from cvx.risk.rand import rand_cov 

62 >>> cov1 = rand_cov(3, seed=42) 

63 >>> cov2 = rand_cov(3, seed=42) 

64 >>> np.allclose(cov1, cov2) 

65 True 

66 

67 Use in portfolio optimization: 

68 

69 >>> from cvx.risk.sample import SampleCovariance 

70 >>> import cvxpy as cp 

71 >>> model = SampleCovariance(num=4) 

72 >>> cov = rand_cov(4, seed=42) 

73 >>> model.update( 

74 ... cov=cov, 

75 ... lower_assets=np.zeros(4), 

76 ... upper_assets=np.ones(4) 

77 ... ) 

78 >>> weights = cp.Variable(4) 

79 >>> risk = model.estimate(weights) 

80 >>> isinstance(risk, cp.Expression) 

81 True 

82 

83 Verify positive definiteness via Cholesky decomposition: 

84 

85 >>> cov = rand_cov(5, seed=123) 

86 >>> # If Cholesky succeeds without error, matrix is positive definite 

87 >>> L = np.linalg.cholesky(cov) 

88 >>> bool(np.allclose(L @ L.T, cov)) 

89 True 

90 

91 Eigenvalue verification: 

92 

93 >>> cov = rand_cov(3, seed=99) 

94 >>> eigenvalues = np.linalg.eigvalsh(cov) 

95 >>> # All eigenvalues should be positive for PD matrix 

96 >>> bool(np.all(eigenvalues > 0)) 

97 True 

98 

99 Different seeds produce different matrices: 

100 

101 >>> cov1 = rand_cov(3, seed=1) 

102 >>> cov2 = rand_cov(3, seed=2) 

103 >>> bool(not np.allclose(cov1, cov2)) 

104 True 

105 

106 Without seed, consecutive calls may differ (random state): 

107 

108 >>> # These may or may not be equal depending on random state 

109 >>> cov_a = rand_cov(2, seed=None) 

110 >>> cov_b = rand_cov(2, seed=None) 

111 >>> cov_a.shape == cov_b.shape == (2, 2) 

112 True 

113 

114 Monte Carlo simulation example: 

115 

116 >>> from cvx.risk.portfolio import minrisk_problem 

117 >>> results = [] 

118 >>> for i in range(5): 

119 ... cov = rand_cov(3, seed=i) 

120 ... model = SampleCovariance(num=3) 

121 ... model.update( 

122 ... cov=cov, 

123 ... lower_assets=np.zeros(3), 

124 ... upper_assets=np.ones(3) 

125 ... ) 

126 ... weights = cp.Variable(3) 

127 ... prob = minrisk_problem(model, weights) 

128 ... _ = prob.solve(solver="CLARABEL") 

129 ... results.append(prob.value) 

130 >>> len(results) 

131 5 

132 >>> all(r > 0 for r in results) # All risks are positive 

133 True 

134 

135 Note: 

136 The generated matrix is guaranteed to be positive semi-definite because 

137 it is constructed as A^T @ A. In practice, it will typically be positive 

138 definite (all eigenvalues strictly positive) unless n is very large. 

139 

140 """ 

141 rng = np.random.default_rng(seed) 

142 a = rng.standard_normal((n, n)) 

143 return np.transpose(a) @ a