Coverage for src/cvx/linalg/decomposition/cholesky.py: 100%

1"""Cholesky decomposition utilities for covariance matrices.""" 

2 

3from __future__ import annotations 

4 

5import warnings 

6from collections.abc import Callable 

7from typing import cast 

8 

9import numpy as np 

10from numpy.linalg import cholesky as _cholesky 

11 

12from ..core.types import Matrix, Vector 

13 

14try: # SciPy's triangular solves keep a factored solve at O(n^2) per right-hand side. 

15 from scipy.linalg import cho_factor as _cho_factor # type: ignore[import-untyped] 

16 from scipy.linalg import cho_solve as _cho_solve 

17 

18 _HAVE_SCIPY = True 

19except ImportError: # pragma: no cover - depends on the environment; the fallback is tested by patching the flag 

20 _HAVE_SCIPY = False 

21 

22 

23def cholesky(cov: Matrix, rhs: Vector | Matrix | None = None) -> Vector | Matrix: 

24 """Compute the upper triangular Cholesky factor of a covariance matrix. 

25 

26 Returns the upper triangular factor R such that R.T @ R = cov. 

27 

28 Args: 

29 cov: A positive definite covariance matrix of shape (n, n). 

30 rhs: Deprecated. When provided the system ``cov @ x = rhs`` is solved 

31 and *x* is returned; use :func:`cholesky_solve` instead. This 

32 parameter will be removed in 2.0. 

33 

34 Returns: 

35 The upper triangular Cholesky factor R when *rhs* is ``None``, or the 

36 solution x to ``cov @ x = rhs`` otherwise (deprecated). 

37 

38 Raises: 

39 np.linalg.LinAlgError: When *rhs* is ``None`` and *cov* is not 

40 positive-definite, or when *rhs* is given and both Cholesky and 

41 LU-based solves fail. 

42 

43 Warns: 

44 DeprecationWarning: When *rhs* is given. 

45 

46 Example: 

47 >>> import numpy as np 

48 >>> from cvx.linalg import cholesky 

49 >>> cov = np.array([[4.0, 2.0], [2.0, 5.0]]) 

50 >>> R = cholesky(cov) 

51 >>> np.allclose(R.T @ R, cov) 

52 True 

53 """ 

54 if rhs is None: 

55 return cast("Matrix", _cholesky(cov).transpose()) 

56 warnings.warn( 

57 "Passing 'rhs' to cholesky() is deprecated and will be removed in 2.0; use cholesky_solve(cov, rhs) instead.", 

58 DeprecationWarning, 

59 stacklevel=2, 

60 ) 

61 return cholesky_solve(cov, rhs) 

62 

63 

64def cholesky_solve(cov: Matrix, rhs: Vector | Matrix) -> Vector | Matrix: 

65 """Solve ``cov @ x = rhs`` using the Cholesky decomposition. 

66 

67 The Cholesky factorisation is attempted first for numerical stability; 

68 when *cov* is not positive-definite the solve falls back to LU 

69 decomposition. With SciPy installed (the ``scipy`` extra) the factor is 

70 applied by two triangular solves, so the whole solve costs one 

71 factorisation. Without it, NumPy has no triangular solve, and a successful 

72 Cholesky is followed by a single LU solve. 

73 

74 Args: 

75 cov: A positive definite covariance matrix of shape (n, n). 

76 rhs: Right-hand side vector of length n or matrix of shape (n, k). 

77 

78 Returns: 

79 The solution x to ``cov @ x = rhs`` with the same shape as *rhs*. 

80 

81 Raises: 

82 np.linalg.LinAlgError: When both the Cholesky and LU-based solves fail. 

83 

84 Example: 

85 >>> import numpy as np 

86 >>> from cvx.linalg import cholesky_solve 

87 >>> cholesky_solve(np.eye(2), np.array([1.0, 2.0])).tolist() 

88 [1.0, 2.0] 

89 >>> cholesky_solve(np.array([[4.0, 0.0], [0.0, 9.0]]), np.array([8.0, 27.0])).tolist() 

90 [2.0, 3.0] 

91 """ 

92 return _factored_solver(cov)(rhs) 

93 

94 

95def _factored_solver(cov: Matrix) -> Callable[[Vector | Matrix], Vector | Matrix]: 

96 """Factorise *cov* once and return a function solving ``cov @ x = rhs``. 

97 

98 The returned solver follows :func:`cholesky_solve` -- Cholesky when *cov* is 

99 positive-definite, LU otherwise -- but pays for the factorisation a single 

100 time, so repeated solves against a fixed matrix cost ``O(n**2)`` each with 

101 SciPy installed. Without SciPy each call is still an LU solve. A matrix that 

102 is not positive-definite falls back to LU at call time, so a singular *cov* 

103 raises :class:`numpy.linalg.LinAlgError` from the solver, not from here. 

104 

105 Args: 

106 cov: A covariance matrix of shape (n, n). 

107 

108 Returns: 

109 A function mapping a right-hand side of shape (n,) or (n, k) to the 

110 solution of the same shape. 

111 """ 

112 try: 

113 if _HAVE_SCIPY and cov.shape[0]: # older SciPy's cho_solve rejects a 0 x 0 factor 

114 factor = _cho_factor(cov) 

115 # cho_factor has already checked cov; skipping the rhs scan lets NaNs propagate. 

116 return lambda rhs: cast("Vector | Matrix", _cho_solve(factor, rhs, check_finite=False)) 

117 _cholesky(cov) # raises LinAlgError unless cov is positive-definite 

118 except np.linalg.LinAlgError: 

119 pass 

120 return lambda rhs: cast("Vector | Matrix", np.linalg.solve(cov, rhs)) 

121 

122 

123def is_positive_definite(matrix: Matrix) -> bool: 

124 """Return True if *matrix* is symmetric positive-definite, False otherwise. 

125 

126 The check is performed via an attempted Cholesky decomposition — the most 

127 numerically reliable way to test positive-definiteness for symmetric matrices. 

128 

129 This function is side-effect-free: it raises no exceptions and emits no 

130 warnings. It is suitable for use as a guard before passing a matrix to a 

131 linear solver. 

132 

133 Args: 

134 matrix: Square matrix to test. 

135 

136 Returns: 

137 ``True`` if the matrix is positive-definite, ``False`` otherwise. 

138 

139 Example: 

140 >>> import numpy as np 

141 >>> from cvx.linalg import is_positive_definite 

142 >>> is_positive_definite(np.eye(3)) 

143 True 

144 >>> is_positive_definite(np.array([[1.0, 2.0], [2.0, 1.0]])) 

145 False 

146 >>> is_positive_definite(np.array([[1.0, 0.5], [0.5, 1.0]])) 

147 True 

148 """ 

149 try: 

150 cholesky(matrix) 

151 except np.linalg.LinAlgError: 

152 return False 

153 else: 

154 return True