Coverage for src/cvx/linalg/operators/dense.py: 100%
126 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-01 07:23 +0000
1"""Dense-matrix backends: :class:`DenseOperator` and :class:`IncrementalDenseOperator`.
3:class:`DenseOperator` wraps an explicit ``n x n`` matrix and slices it directly.
4:class:`IncrementalDenseOperator` specialises it for active-set sweeps, maintaining
5the free-block inverse across single-index changes with in-place rank-one bordered /
6deletion updates instead of refactorising each step.
7"""
9from __future__ import annotations
11import numpy as np
13from ..core.exceptions import DimensionMismatchError, NonSquareMatrixError, NotAMatrixError
14from ..core.types import Matrix, Vector
15from ..decomposition.cholesky import cholesky_solve
16from .base import SymmetricOperator, as_index, rcond_symmetric
18try: # BLAS ``ger`` applies the rank-one updates in place, with no k x k temporary.
19 from scipy.linalg.blas import get_blas_funcs as _get_blas_funcs # type: ignore[import-untyped]
21 _HAVE_SCIPY = True
22except ImportError: # pragma: no cover - depends on the environment; the fallback is tested by patching the flag
23 _HAVE_SCIPY = False
26class DenseOperator(SymmetricOperator):
27 """Symmetric operator backed by an explicit dense matrix.
29 Args:
30 matrix: A symmetric ``n x n`` matrix. It is stored by reference, not
31 copied or symmetrised.
33 Example:
34 >>> import numpy as np
35 >>> from cvx.linalg import DenseOperator
36 >>> A = np.array([[4.0, 1.0, 0.0], [1.0, 3.0, 1.0], [0.0, 1.0, 2.0]])
37 >>> op = DenseOperator(A)
38 >>> free = np.array([0, 2])
39 >>> np.allclose(op.apply_free(free, np.array([1.0, 1.0])), A[np.ix_(free, free)] @ np.ones(2))
40 True
41 """
43 def __init__(self, matrix: Matrix) -> None:
44 """Store the backing matrix after checking it is square."""
45 matrix = np.asarray(matrix, dtype=np.float64)
46 if matrix.ndim != 2:
47 raise NotAMatrixError(matrix.ndim, func="DenseOperator")
48 if matrix.shape[0] != matrix.shape[1]:
49 raise NonSquareMatrixError(matrix.shape[0], matrix.shape[1])
50 self._a = matrix
52 @property
53 def n(self) -> int:
54 """Dimension of the operator."""
55 return int(self._a.shape[0])
57 @property
58 def diag(self) -> Vector:
59 """The diagonal of the backing matrix (a read-only view)."""
60 result: Vector = np.diagonal(self._a)
61 return result
63 def matvec(self, x: Vector | Matrix) -> Vector | Matrix:
64 """Return ``A @ x`` by dense multiplication."""
65 return self._a @ x
67 def restricted(self, free: object) -> DenseOperator:
68 """Return ``DenseOperator(A[free, free])``: the free block, pre-sliced."""
69 free = as_index(free)
70 return DenseOperator(self._a[np.ix_(free, free)])
72 def block_matvec(self, rows: object, cols: object, v: Vector | Matrix) -> Vector | Matrix:
73 """Return ``A[rows, cols] @ v`` by slicing the dense matrix."""
74 rows = as_index(rows)
75 cols = as_index(cols)
76 result: Vector | Matrix = self._a[np.ix_(rows, cols)] @ v
77 return result
79 def solve_free(self, free: object, rhs: Vector | Matrix) -> Vector | Matrix:
80 """Solve the free block by Cholesky (LU fallback via :func:`cholesky_solve`)."""
81 free = as_index(free)
82 return cholesky_solve(self._a[np.ix_(free, free)], rhs)
84 def rcond_free(self, free: object) -> float:
85 """Reciprocal condition number of the free block from its symmetric eigenvalues."""
86 free = as_index(free)
87 return rcond_symmetric(self._a[np.ix_(free, free)])
90class IncrementalDenseOperator(DenseOperator):
91 """Dense operator that maintains ``A[free, free]^{-1}`` across single-index flips.
93 A drop-in :class:`DenseOperator` whose :meth:`solve_free` reuses the previous
94 free-block inverse when the free set changed by at most one index since the last
95 call, updating it with a rank-one bordered (index added) or deletion (index
96 removed) formula at ``O(n * len(free))`` instead of refactorising at
97 ``O(len(free)**3)``. Any other change -- a multi-index change, or a
98 non-positive/non-finite pivot -- recomputes the inverse from scratch.
100 The inverse is kept in *insertion order* in the leading block of a preallocated
101 ``n x n`` Fortran-ordered buffer (allocated on the first solve), so an update
102 never copies or permutes the whole block: an insert appends a border row and
103 column, and a delete swaps the leaving slot with the last one. With SciPy
104 installed the rank-one term is applied in place by BLAS ``ger``; without it, by
105 NumPy with one ``k x k`` temporary.
107 This suits an active-set loop that changes its free set one index at a time. The
108 free indices may come in any order; *rhs* is aligned to that order, and so is the
109 returned solution. :meth:`matvec`, :meth:`block_matvec`, and :meth:`rcond_free`
110 are the plain dense ones -- only :meth:`solve_free` differs.
112 A maintained inverse accumulates rounding over the ``O(n)`` updates of a sweep, so
113 on ill-conditioned problems the plain :class:`DenseOperator` (a clean solve each
114 step) is the safer choice.
116 Example:
117 >>> import numpy as np
118 >>> from cvx.linalg import IncrementalDenseOperator
119 >>> op = IncrementalDenseOperator(np.eye(3))
120 >>> np.allclose(op.solve_free(np.array([0, 1]), np.array([1.0, 2.0])), [1.0, 2.0])
121 True
122 """
124 def __init__(self, matrix: Matrix) -> None:
125 """Wrap ``matrix`` (validated as in :class:`DenseOperator`) and start with no cache."""
126 super().__init__(matrix)
127 n = self.n
128 self._range = np.arange(n, dtype=np.intp) # bounds-checks and normalises free indices
129 self._order = np.empty(n, dtype=np.intp) # index held in each slot
130 self._slot = np.full(n, -1, dtype=np.intp) # slot of each index, -1 when not free
131 self._k = 0 # number of occupied slots
132 self._buf: Matrix = np.empty((0, 0), order="F") # inverse in slot order in its leading k x k block
134 def solve_free(self, free: object, rhs: Vector | Matrix) -> Vector | Matrix:
135 """Solve ``A[free, free] @ y = rhs`` using the maintained (incrementally updated) inverse."""
136 cur = self._range[as_index(free)]
137 rhs = np.asarray(rhs)
138 if rhs.ndim == 0 or rhs.shape[0] != cur.size:
139 raise DimensionMismatchError(rhs.shape[0] if rhs.ndim else rhs.size, cur.size)
140 if self._buf.shape[0] != self.n: # allocated on the first solve
141 self._buf = np.zeros((self.n, self.n), order="F")
142 if not self._update(cur):
143 self._refactor(cur)
144 pos = self._slot[cur]
145 rhs_slot = np.empty(rhs.shape, dtype=np.result_type(rhs, self._buf))
146 rhs_slot[pos] = rhs
147 solution: Vector | Matrix = (self._buf[: self._k, : self._k] @ rhs_slot)[pos]
148 return solution
150 def _update(self, cur: np.ndarray) -> bool:
151 """Bring the cache from the held free set to *cur* by at most one update (``False`` if it cannot)."""
152 n_held = int(np.count_nonzero(self._slot[cur] >= 0))
153 n_added, n_removed = cur.size - n_held, self._k - n_held
154 if n_added == 0 and n_removed == 0:
155 return True
156 if n_added == 1 and n_removed == 0:
157 return self._insert(int(cur[self._slot[cur] < 0][0]))
158 if n_added == 0 and n_removed == 1:
159 held = self._order[: self._k]
160 keep = np.zeros(self.n, dtype=bool)
161 keep[cur] = True
162 return self._delete(int(held[~keep[held]][0]))
163 return False # not a single-index flip; recompute
165 def _refactor(self, cur: np.ndarray) -> None:
166 """Invert ``A[cur, cur]`` from scratch into slots ``0..len(cur)-1``, in *cur* order."""
167 k = cur.size
168 inv = np.linalg.inv(self._a[np.ix_(cur, cur)]) # may raise: leave the cache intact
169 self._slot[self._order[: self._k]] = -1
170 self._order[:k] = cur
171 self._slot[cur] = self._range[:k]
172 self._k = k
173 self._block()[:] = inv
175 def _block(self) -> Matrix:
176 """The leading ``k x k`` block of the buffer: the maintained inverse in slot order."""
177 block: Matrix = self._buf[: self._k, : self._k]
178 return block
180 def _rank_one(self, alpha: float, x: Vector) -> None:
181 """Add ``alpha * outer(x, x)`` to the leading ``len(x) x len(x)`` block of the buffer, in place."""
182 m = x.shape[0]
183 buf = self._buf
184 if m == 0:
185 return
186 if _HAVE_SCIPY:
187 # buf[:, :m] is Fortran-contiguous, so ger updates it in place; buf[:m, :m] is
188 # not, and ger would copy it. Zero-padding x leaves rows m.. unchanged.
189 padded = np.zeros(self.n)
190 padded[:m] = x
191 ger = _get_blas_funcs("ger", (buf,))
192 ger(alpha, padded, x, a=buf[:, :m], overwrite_a=True)
193 else:
194 block = buf[:m, :m]
195 block += alpha * np.outer(x, x)
197 def _insert(self, asset: int) -> bool:
198 """Rank-one bordered update for one index entering the free set (``False`` if the pivot is bad)."""
199 k = self._k
200 c = self._a[self._order[:k], asset]
201 v = self._block() @ c
202 schur = float(self._a[asset, asset] - c @ v)
203 if not np.isfinite(schur) or schur <= 0.0:
204 return False
205 self._rank_one(1.0 / schur, v)
206 buf = self._buf
207 buf[:k, k] = buf[k, :k] = -v / schur
208 buf[k, k] = 1.0 / schur
209 self._order[k], self._slot[asset], self._k = asset, k, k + 1
210 return True
212 def _delete(self, asset: int) -> bool:
213 """Rank-one deletion update for one index leaving the free set (``False`` if the pivot is bad)."""
214 block = self._block()
215 p, last = int(self._slot[asset]), self._k - 1
216 pivot = float(block[p, p])
217 if not np.isfinite(pivot) or pivot <= 0.0:
218 return False
219 if p != last: # move the index in the last slot into slot p
220 block[[p, last], :] = block[[last, p], :]
221 block[:, [p, last]] = block[:, [last, p]]
222 moved = self._order[last]
223 self._order[p], self._slot[moved] = moved, p
224 self._rank_one(-1.0 / pivot, block[:last, last].copy())
225 self._slot[asset], self._k = -1, last
226 return True