Coverage for src/cvx/linalg/operators/base.py: 100%

1"""The :class:`SymmetricOperator` protocol and shared index/conditioning helpers. 

2 

3This module holds the abstract contract every backend implements together with 

4the two structure-agnostic helpers the backends share: :func:`_as_index`, which 

5normalises an index set, and :func:`_rcond_symmetric`, which reads the reciprocal 

62-norm condition number off a symmetric block's eigenvalues. 

7""" 

8 

9from __future__ import annotations 

10 

11from abc import ABC, abstractmethod 

12 

13import numpy as np 

14 

15from ..core.types import Matrix, Vector 

16 

17_RESTRICTED_DEFAULT_MESSAGE = "this backend does not implement restricted(free)" 

18_INDEX_NDIM_MESSAGE = "index set must be a 1-D array of integer positions" 

19_INDEX_DUPLICATE_MESSAGE = "index set must not contain duplicates" 

20_NO_DIAG_MESSAGE = "this operator does not expose its diagonal" 

21 

22 

23def _as_index(indices: object) -> np.ndarray: 

24 """Coerce an index set to a 1-D integer array.""" 

25 arr = np.asarray(indices) 

26 if arr.ndim != 1: 

27 raise ValueError(_INDEX_NDIM_MESSAGE) 

28 if arr.dtype == np.bool_ or not np.issubdtype(arr.dtype, np.integer): 

29 raise ValueError(_INDEX_NDIM_MESSAGE) 

30 idx = arr.astype(np.intp, copy=False) 

31 if np.unique(idx).size != idx.size: 

32 raise ValueError(_INDEX_DUPLICATE_MESSAGE) 

33 return idx 

34 

35 

36def _rcond_symmetric(block: Matrix) -> float: 

37 """Reciprocal 2-norm condition number of a symmetric positive-(semi)definite block. 

38 

39 Returns a value in ``[0, 1]`` read off the symmetric eigenvalues 

40 (``lambda_min / lambda_max``): ``1`` is perfectly conditioned, ``0`` numerically 

41 singular. Deterministic and BLAS-independent. An empty block is well conditioned. 

42 """ 

43 if block.shape[0] == 0: 

44 return 1.0 

45 eigenvalues = np.linalg.eigvalsh(block) 

46 lam_max = float(eigenvalues[-1]) 

47 if lam_max <= 0.0: 

48 return 0.0 

49 return max(float(eigenvalues[0]), 0.0) / lam_max 

50 

51 

52class SymmetricOperator(ABC): 

53 """A symmetric operator accessed only through block products and a free solve. 

54 

55 Subclasses implement :meth:`matvec`, :meth:`block_matvec`, and 

56 :meth:`solve_free`; :meth:`apply_free` is derived and :attr:`diag` is 

57 optional (it raises unless the backend overrides it). Index sets passed to 

58 the block methods are 1-D arrays of integer positions into ``range(n)``. 

59 """ 

60 

61 @property 

62 @abstractmethod 

63 def n(self) -> int: 

64 """Dimension of the operator (it acts on vectors of length ``n``).""" 

65 

66 @property 

67 def diag(self) -> Vector: 

68 """The diagonal of ``A`` as a length-``n`` vector. 

69 

70 The natural Jacobi preconditioner for a Krylov solve over 

71 :meth:`apply_free` (the free block's diagonal is ``diag[free]``). 

72 Backends override this with a closed form read off their structure; 

73 the default raises for operators without a cheap diagonal. 

74 

75 Raises: 

76 NotImplementedError: If the backend does not expose its diagonal. 

77 """ 

78 raise NotImplementedError(_NO_DIAG_MESSAGE) 

79 

80 @abstractmethod 

81 def matvec(self, x: Vector | Matrix) -> Vector | Matrix: 

82 """Return ``A @ x``.""" 

83 

84 @abstractmethod 

85 def block_matvec(self, rows: object, cols: object, v: Vector | Matrix) -> Vector | Matrix: 

86 """Return ``A[rows, cols] @ v``, the product of a sub-block with *v*. 

87 

88 Here ``A[rows, cols]`` is the sub-matrix with row indices *rows* and 

89 column indices *cols* (as by ``A[np.ix_(rows, cols)]``); *v* has length 

90 ``len(cols)`` (or that many rows) and the result has length 

91 ``len(rows)``. 

92 """ 

93 

94 @abstractmethod 

95 def solve_free(self, free: object, rhs: Vector | Matrix) -> Vector | Matrix: 

96 """Return ``A[free, free]^{-1} @ rhs``, a direct solve on a principal block. 

97 

98 ``A[free, free]`` must be positive definite; a Krylov solver that would 

99 rather iterate should call :meth:`apply_free` instead. 

100 """ 

101 

102 @abstractmethod 

103 def rcond_free(self, free: object) -> float: 

104 """Return the reciprocal 2-norm condition number of ``A[free, free]``. 

105 

106 A value in ``[0, 1]`` (``1`` well conditioned, ``0`` numerically singular). 

107 An active-set solver uses this to detect a rank-deficient free block before 

108 trusting :meth:`solve_free`. Backends compute it from their structure 

109 without materialising the ``len(free) x len(free)`` block. 

110 """ 

111 

112 def restricted(self, free: object) -> SymmetricOperator: 

113 """Return the operator of ``A[free, free]`` with pre-sliced storage. 

114 

115 Use this to hoist the free-set restriction out of an iterative solve: 

116 build the restricted operator once per working set, then call its 

117 :meth:`matvec` inside the loop. Calling :meth:`apply_free` per 

118 iteration instead re-gathers the underlying storage (e.g. the factor 

119 columns of a Gram operator) on every call, which can dominate the 

120 useful arithmetic by an order of magnitude. Backends override this 

121 with a pre-sliced copy of their structure; the default raises. 

122 """ 

123 raise NotImplementedError(_RESTRICTED_DEFAULT_MESSAGE) 

124 

125 def apply_free(self, free: object, v: Vector | Matrix) -> Vector | Matrix: 

126 """Return ``A[free, free] @ v``, the free-block forward product. 

127 

128 It is :meth:`block_matvec` with equal row and column sets. Each call 

129 gathers the free-block storage; for repeated products on a fixed free 

130 set (a CG loop), build :meth:`restricted` once and use its 

131 :meth:`matvec`. 

132 """ 

133 free = _as_index(free) 

134 return self.block_matvec(free, free, v)