Coverage for src/nncg/certificate.py: 100%

1"""The KKT certificate for the non-negative quadratic program and its shared precondition. 

2 

3:func:`kkt_violation` scores how far a candidate is from the unique global 

4minimiser of ``min_{x>=0} 1/2 x'Ax - b'x`` — zero certifies optimality — and is 

5the load-bearing check the paper's numerical study reports against. 

6:func:`_require_operator` is the one operator/right-hand-side precondition shared 

7by the certificate and both :class:`nncg.solver.ActiveSetSolver` entry points. 

8""" 

9 

10from __future__ import annotations 

11 

12import numpy as np 

13from cvx.linalg import SymmetricOperator, Vector 

14 

15_NEEDS_OPERATOR = ( 

16 "the quadratic term must be a cvx.linalg.SymmetricOperator: wrap a dense SPD " 

17 "array in DenseOperator(a), or pass GramOperator(M, ridge) for A = M'M + ridge*I" 

18) 

19 

20 

21def _require_operator(a: SymmetricOperator, b: Vector) -> None: 

22 """Validate that ``a`` is a symmetric operator whose dimension matches ``b``. 

23 

24 Args: 

25 a: The quadratic term, expected to be a :class:`cvx.linalg.SymmetricOperator`. 

26 b: The linear term ``b``. 

27 

28 Raises: 

29 TypeError: When ``a`` is not a :class:`cvx.linalg.SymmetricOperator`. 

30 ValueError: When the operator dimension does not match ``len(b)``. 

31 """ 

32 if not isinstance(a, SymmetricOperator): 

33 raise TypeError(_NEEDS_OPERATOR) 

34 if a.n != len(b): 

35 msg = f"operator dimension {a.n} does not match len(b) = {len(b)}" 

36 raise ValueError(msg) 

37 

38 

39def kkt_violation(a: SymmetricOperator, b: Vector, x: Vector) -> float: 

40 """Maximum violation of the KKT system of ``min_{x>=0} 1/2 x'Ax - b'x``. 

41 

42 Args: 

43 a: The SPD operator ``A`` (a :class:`cvx.linalg.SymmetricOperator`). 

44 b: The linear term ``b``. 

45 x: Candidate solution. 

46 

47 Returns: 

48 ``max`` of the negativity violations of ``x`` and of the reduced 

49 gradient ``s = A x - b``, and of the complementarity products 

50 ``|x_i s_i|``. Zero certifies the unique global minimiser. 

51 

52 Examples: 

53 Note that ``a`` must be an operator — a bare array raises ``TypeError``: 

54 

55 >>> import numpy as np 

56 >>> from cvx.linalg import DenseOperator 

57 >>> a = DenseOperator(np.array([[2.0, 0.0], [0.0, 2.0]])) 

58 >>> b = np.array([2.0, -2.0]) 

59 

60 The minimiser certifies at zero, while the origin does not: 

61 

62 >>> round(kkt_violation(a, b, np.array([1.0, 0.0])), 12) 

63 0.0 

64 >>> kkt_violation(a, b, np.zeros(2)) > 0 

65 True 

66 """ 

67 _require_operator(a, b) 

68 s = a.matvec(x) - b 

69 # The leading 0.0 is a floor on a quantity that is mathematically non-negative 

70 # already, and it is there for the sign of zero rather than the magnitude: when 

71 # every term is zero, `np.max(-s, initial=0.0)` may hand back -0.0 (its reduce 

72 # path picks a different one of the two equal zeros on linux than on macOS), and 

73 # `-0.0` compares equal to `0.0` but does not print the same. Python's `max` 

74 # replaces its running best only on a strict `>`, so the +0.0 seeded here 

75 # survives and the certificate reports one platform-independent zero. 

76 return float( 

77 max( 

78 0.0, 

79 np.max(-x, initial=0.0), 

80 np.max(-s, initial=0.0), 

81 np.max(np.abs(x * s), initial=0.0), 

82 ) 

83 )