Coverage for src/cvx/quadprog/_setup.py: 100%

1"""Everything that happens before the first iteration. 

2 

3Validation of the caller's arrays, and the factorisation of ``G`` that the 

4iteration is carried out in terms of. 

5""" 

6 

7# G, C, R and J are the names used in Goldfarb & Idnani (1983) and in the 

8# reference implementation's public signature `solve_qp(G, a, C, b, meq)`. 

9# Lowercasing them would obscure the correspondence to the paper, so the 

10# pep8-naming rules are waived here, as they are in _solve.py. 

11# ruff: noqa: N803, N806, TRY003 

12 

13import numpy as np 

14import scipy.linalg 

15from scipy.linalg.lapack import dtrtri 

16 

17 

18def _validate( 

19 G: np.ndarray, a: np.ndarray, C: np.ndarray, b: np.ndarray, meq: int, check_finite: bool = False 

20) -> tuple[int, int]: 

21 """Check that the problem data is dimensionally consistent. 

22 

23 Args: 

24 G: ``(n, n)`` matrix of the quadratic term. 

25 a: ``(n,)`` vector of the linear term. 

26 C: ``(n, m)`` constraint matrix. 

27 b: ``(m,)`` right-hand side of the constraints. 

28 meq: Number of leading constraints treated as equalities. 

29 check_finite: Whether to reject NaN and infinity in the inputs. Off by 

30 default, matching the reference; see ``solve_qp``. 

31 

32 Returns: 

33 The number of variables and the number of constraints. 

34 

35 Raises: 

36 ValueError: If any shape disagrees, if ``meq`` is out of range, or if 

37 ``check_finite`` is set and any input holds a non-finite value. 

38 """ 

39 if G.ndim != 2 or G.shape[0] != G.shape[1]: 

40 raise ValueError(f"G must be a square matrix. Received shape={G.shape}") 

41 n = G.shape[0] 

42 if a.shape != (n,): 

43 raise ValueError(f"G and a must have the same dimension. Received G as {G.shape} and a as {a.shape}") 

44 if C.ndim != 2 or C.shape[0] != n: 

45 raise ValueError(f"G and C must have the same first dimension. Received G as {G.shape} and C as {C.shape}") 

46 q = C.shape[1] 

47 if b.shape != (q,): 

48 raise ValueError( 

49 f"The number of columns of C must match the length of b. Received C as {C.shape} and b as {b.shape}" 

50 ) 

51 if not 0 <= meq <= q: 

52 raise ValueError(f"meq must satisfy 0 <= meq <= {q}. Received {meq}") 

53 if check_finite: 

54 # Last, so a caller who passes both a wrong shape and a NaN still hears 

55 # about the shape -- that is the error they can act on without reading 

56 # their data. 

57 _check_finite(G, a, C, b) 

58 return n, q 

59 

60 

61def _check_finite(G: np.ndarray, a: np.ndarray, C: np.ndarray, b: np.ndarray) -> None: 

62 """Reject NaN and infinity in the problem data, naming the first offender. 

63 

64 Only reached when ``check_finite`` is set: the scan is :math:`O(n^2)` on 

65 ``G``, which is why it is opt-in rather than unconditional. 

66 

67 Args: 

68 G: ``(n, n)`` matrix of the quadratic term. 

69 a: ``(n,)`` vector of the linear term. 

70 C: ``(n, m)`` constraint matrix. 

71 b: ``(m,)`` right-hand side of the constraints. 

72 

73 Raises: 

74 ValueError: If any argument holds a non-finite value. 

75 """ 

76 for name, array in (("G", G), ("a", a), ("C", C), ("b", b)): 

77 if not np.isfinite(array).all(): 

78 raise ValueError(f"{name} contains a non-finite value (NaN or infinity)") 

79 

80 

81def _factorize(G: np.ndarray, a: np.ndarray, factorized: bool) -> tuple[np.ndarray, np.ndarray]: 

82 """Return the inverse Cholesky factor of ``G`` and the unconstrained minimum. 

83 

84 Args: 

85 G: ``(n, n)`` positive definite matrix, or its inverse Cholesky factor 

86 :math:`R^{-1}` when ``factorized`` is True. 

87 a: ``(n,)`` vector of the linear term. 

88 factorized: Whether ``G`` already holds :math:`R^{-1}`. 

89 

90 Returns: 

91 ``J``, an upper triangular array with ``J J^T = G^-1``, and the 

92 unconstrained minimiser ``G^-1 a``. ``J`` is a fresh writable array; the 

93 caller updates it in place. 

94 

95 Raises: 

96 ValueError: If ``G`` is not positive definite. 

97 """ 

98 # Fortran order throughout: the updates in _qr work on column blocks of J, 

99 # which are then contiguous and can go straight to BLAS. 

100 if factorized: 

101 J = np.asfortranarray(np.triu(G)) 

102 return J, J @ (J.T @ a) 

103 

104 # check_finite=False skips a full scan of each array on the way in, which 

105 # matches the reference: it does not check either. A non-finite G is 

106 # therefore not diagnosed here, and what happens next is a property of the 

107 # LAPACK build rather than of this package -- Accelerate reports a failed 

108 # potrf and raises below, OpenBLAS runs to completion and propagates NaNs 

109 # into the result. Neither returns a finite wrong answer, which is the only 

110 # guarantee callers can portably rely on. 

111 try: 

112 R = scipy.linalg.cholesky(G, lower=False, check_finite=False) 

113 except scipy.linalg.LinAlgError as exc: 

114 raise ValueError("matrix G is not positive definite") from exc 

115 

116 xv = scipy.linalg.cho_solve((R, False), a, check_finite=False) 

117 J, info = dtrtri(R, lower=0) 

118 if info != 0: # pragma: no cover 

119 # Defensive: trtri fails only on an exactly zero diagonal entry, which 

120 # a successful Cholesky has already ruled out. 

121 raise ValueError("matrix G is not positive definite") 

122 return np.asfortranarray(np.triu(J)), xv