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

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

2 

3from __future__ import annotations 

4 

5import numpy as np 

6from numpy.linalg import cholesky as _cholesky 

7 

8 

9def cholesky(cov: np.ndarray, rhs: np.ndarray | None = None) -> np.ndarray: 

10 """Compute the upper triangular Cholesky factor, or solve a linear system. 

11 

12 When called without *rhs*, returns the upper triangular factor R such that 

13 R.T @ R = cov. When *rhs* is given, solves ``cov @ x = rhs`` using the 

14 Cholesky decomposition for numerical stability, falling back to LU 

15 decomposition if *cov* is not positive-definite. 

16 

17 Args: 

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

19 rhs: Optional right-hand side vector or matrix. When provided the 

20 system ``cov @ x = rhs`` is solved and *x* is returned. 

21 

22 Returns: 

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

24 solution x to ``cov @ x = rhs`` otherwise. 

25 

26 Raises: 

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

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

29 LU-based solves fail. 

30 

31 Example: 

32 Decomposition (no rhs): 

33 

34 >>> import numpy as np 

35 >>> from cvx.linalg import cholesky 

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

37 >>> R = cholesky(cov) 

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

39 True 

40 

41 Solve (with rhs): 

42 

43 >>> cholesky(np.eye(2), np.array([1.0, 2.0])).tolist() 

44 [1.0, 2.0] 

45 >>> cholesky(np.array([[4.0, 0.0], [0.0, 9.0]]), np.array([8.0, 27.0])).tolist() 

46 [2.0, 3.0] 

47 """ 

48 if rhs is None: 

49 return _cholesky(cov).transpose() 

50 try: 

51 upper = _cholesky(cov).transpose() 

52 return np.linalg.solve(upper, np.linalg.solve(upper.T, rhs)) 

53 except np.linalg.LinAlgError: 

54 return np.linalg.solve(cov, rhs) 

55 

56 

57def is_positive_definite(matrix: np.ndarray) -> bool: 

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

59 

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

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

62 

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

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

65 linear solver. 

66 

67 Args: 

68 matrix: Square matrix to test. 

69 

70 Returns: 

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

72 

73 Example: 

74 >>> import numpy as np 

75 >>> from cvx.linalg import is_positive_definite 

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

77 True 

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

79 False 

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

81 True 

82 """ 

83 try: 

84 cholesky(matrix) 

85 except np.linalg.LinAlgError: 

86 return False 

87 else: 

88 return True